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I  .   INTRODUCTION 

A.   HISTORY 

Since  man  first  poured  a  colored  liquid  into  water  and 
observed  that  the  liquid  spread  throughout  the  container, 
since  he  placed  two  metals  together  in  heat  and  later 
observed  that  they  melded  together,  or  since  he  first 
noticed  that  the  smoke  from  his  fire  mixed  into  the  air, 
diffusion  has  been  observed  and  its  cause  a  source  of 
wond  er . 

Initial  qualitative  experiments  in  diffusion  were 
conducted  by  Parrot  in  1815  (Crank,  1975).  In  this 
experiment,  he  noted  that  no  matter  how  carefully  he 
controlled  mechanical  agitation  and  convection,  the  separate 
gases  in  the  experiment  always  diffused  throughout  the 
container.  One  other  important  precursor  development  was 
Fourier's  development  of  the  equation  for  heat  conduction 
(Fourier,  13  22). 

The  history  of  the  scientific  study  of  diffusion  is 
rooted  in  a  paper  by  Adolf  Fick  (IP  5  5  ) .  F  i  c  k  used  Fourier's 
mathematical  methods  to  describe  the  equations  of 
diffusion.  One  equation  states  that  the  rate  of  diffusion 
is  proportional  to  the  concentration  gradient  across  the 
plane  of  diffusion: 

J  =-  D(3  C/  dx),  (  1) 
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w  h  ere  J  is  t  h  e  f  1  u  x  of  diffusing  materia].,  i)  t  h  e  d  i  f  f  u  s  i  o  n 
coefficient,  and  C  the  -concentration  of  the  d  i  f  f  u  s  i 
material.  The  units  of  the  diffusion  coefficient  are 
(length)  /(time),  and  it  is  usually  given  in  centimeters 
per  second.  Equation  (1)  is  known  as  Fick's  first  law  of 
diffusion.  In  this  work,  Equation  (1),  is  assumed  to  be  of 
this  form,  even  for  non-constant  D . 

Fick  also  developed  a  second  equation  which  describes 
the  time  rate  of  change  of  the  concentration  of  the 
diffusing  material.  Known  as  Fick's  second  law,  it  states 
that  the  partial  derivative  of  the  concentration  with 
respect  to  time  is  equal  to  the  negative  of  the  divergence 
of  the  flux: 

3C/3t  =  -V  J.  (2) 

In  one  dimensional  cartesion  coordinates,  this  reduces  to: 


3C/  3t  =  -  3J/  3x. 


(3) 


When  Equation  (3)  and  Fick's  first  law  (Equation  1)  are 
combined  Equation  (3)  becomes: 


3C/  3t  =  {3  /  3x[  (D)(  3C/  3x)  ]  }  , 


(3a) 


and,  if  one  assumes  that  the  diffusion  coefficient  is 
constant  with  respect  to  position,  the  equation  further 
reduces  to: 


3C/3  t  =  D(  32C/  3x2)  . 


(3b) 


B 


The  general  solution  to  this  equation,  given  an  initial 
impulsive  planar  source  which  then  spreads  across  a  one- 
dimensional  infinite  material,  is: 


C  =  (A/t1/2)exp(-x2/4Dt) , 


(A) 


where  A  is  a  constant,  t  the  elapsed  time,  and  x  the 
position  of  interest  within  the  material.  A  consequence  of 
this  solution  is ' that  after  sufficient  time  the 
concentration  will  become  uniform  throughout  the  material  as 
seen  in  Figure  1  . 


Figure  1 .   Concentration  Profiles  at  Different  Times 


Einstein  (1905)  provided  an  atomistic  definition  of  the 
diffusion  coefficient: 

D  =  <x2>/2t  ,  (5) 

where  _x  is  the  vector  displacement  of  a  given  atom  from  its 
initial  position  to  its  final  position  at  time  t.  The  Dim c 
b rackets  denote  an  average  over  a  large  number  of  diffusing 


atoms.   It  can  readily  be  s  h  o  w  n  that  the  -definitions  for  D 

in  Equations  (1)  and  (5)  are  equivalent. 

B.   DIFFUSION  MECHANISMS 

Since  Fick's  and  Einstein's  work  in  diffusion  theory, 
many  people  have  tried  to  explain  more  fully  trie  mechanisms 
by  which  an  atom  diffuses  through  another  material.  The 
majority  of  this  work  was  done  by  Kirkendall  (1942),  H  a  r  1 1 y 
(19  4  6)  and  Darken  (1948).  These  men  developed  tools  to 
explain  diffusion  based  on  a  knowledge  of  the  nature  of  the 
atom  and  the  method  by  which  atoms  combine  into  bulk 
materials.  Given  the  concept  of  the  crystal  structure?  of 
solids,  many  mechanisms  of  diffusion  can  be  identified. 
Table  1  is  a  partial  list  of  the  mechanisms  of  diffusion 
through  solid  crystals: 

TABLE  1 
PARTIAL  TABLE  OF  DIFFUSION  MECHANISMS 

1 .  Interstitial  mechanism 

2.  Inter  stitialcy  mechanism 

3.  Crowdion  mechanism 

4 .  Vacancy  mechanism 

5.  Relax  ion  mechanism. 

The  INTERSTITIAL  MECHANISM  describes  the  diffusion  of 
small  atoms  through  a  crystal  structure  of  larger  atoms.  The 
a  t  o m  advances  thro  u g h  the  structure  b y  m  o v  i n  g  from  o  n  e 
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vacant  Lntertitial  site  within  the  crystal  to  another.   This 
is  shown  in  Figure  2 . 


Figure  2.   lntertitial  Mechanism  of  Diffusion. 

When  the  interstitial  atom's  size  approaches  that  of  the 
lattice  atoms,  they  move  by  the  INTEP.STITI ALCY  MECHANISM.  In 
this  method  the  defect  atom  moves  into  a  normal  lattice 
position  by  pushing  the  lattice  atom  out  into  an 
interstitial  site.   Figure  3  illustrates  this  mechanism. 


Figure  3.   Intertitialcy  Mechanism  of  Diffusion. 


The  C R  0  W DIP  N  MECHANISM  of  diffusion  requires  a  defect 
atom  with  an  energy  large  enough  such  that  it  can  push  L  t  s 
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way  onto  a  lattice  site, compressing  a  close  packed  row  of 
atoms.    This  mechanism  can  be  thought  of  as  a   line 
billiard  balls  touching  each  other.   Then  a  cue  ball  hits 
this  line  so  that  each  ball  moves  one  position,  replacing 
its  neighbor.   An  example  is  shown  in  Figure  4 . 


Figure  4.   Crowd  ion  Mechanism  of  Diffusion. 

All  crystals  have  vacancies  in  their  lattices  at  a 
finite  temperature.  In  the  VACANCY  MECHANISM  of  diffusion, 
these  vacancies  are  filled  by  the  diffusing  atoms.  This 
method  of  diffusion  is  predominant  in  self  diffusion  of  a 
material  since  a  lattice  atom  and  the  vacancy  simply 
exchange  position.  Figure  5  illustrates  the  method  by  which 
the  vacancy  moves. 


Figure  5.   Vacancy  M echanism  of  Diffusion 
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The  RELAX  I ON   MECHANISM  of  diffusion   is  a   modified 


vacancy  mechani-sn 


In  the  area  of  a  vacancy  t  h  ^  lattice 


atoms  relax  towards  the  vacancy.  These  atoms,  and  any  defect 
atoms  in  this  area,  are  able  to  nove  by  an  irregular 
movement  similar  to  motion  of  a  liquid.  This  mechanism  can 
be  seen  in  Figure  6. 


Figure  6.   Relaxion  Mechanism  of  Diffusion. 

J.R.  Manning  (1968)  has   an  excellent  discussion  of  the 
mechanisms  of  diffusion. 

C.   BASIC  CONCEPTS 

Crystaline  solids  are  made  up  of  atoms  that  arrange 
themselves  in  various  lattice  structures.  There  are  seven 
types  of  crystals,   these  are: 


TABLE  2 

TABLE  Or  CRYSTAL  STRUCTURE 

1  .  cubic 

2.  tetragonal 

3.  orthorhonbic 
4  .  hexagona 1 

5.  rhombo'neclral 

6  .  nionoclinic 

7.  triclinic 

Within  each  of  these  types  of  crystals  subdivisions  exist 
that  further  divide  these  seven  structures  into  a  total  of 
fourteen  crystal  structures.  L.A.  Girifalco  (19 6 A)  gives  a 
good  discussion  of  the  different  crystal  structure  s .  The 
face-centered  cubic  (fee)  crystal  structure  is  used  in  this 
thesis  . 


(1  0  0)  (3-D) 

Figure  7 .   Different  Propectives  of  the  fee  Crystal. 

Figure  7  illustrates  different  sections  of  a  fee  crystal. 
There  exists  sites  between  the  a  tons  of  the  structure  known 
as  interstitial  sites,  see  Figure  8 . 
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Figure    8 .    Intertitial    and    Lattice    Sites     of     fee     Crystal. 

Each  site  in  the  fee  structure,  whether  an  interstitial  site 
or  lattice  site,  has  twelve  nearest  neighbors,  Figure  9 
shows  the  relationship  of  a  site  and  its  nearest  neighbors 
within  the  crystal.  The  interstitial  method  of  diffusion 
described  above  moves  an  atom  to  one  of  these  nearest- 
neiehbor     sites. 


(2    -    interstitial     sites 


-    lattice    sites 


i g u r e    9 .       Fee    Atom    and     its    Twelve    Clearest    Neighbors 
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The    force    acting    bet w e e n    a     pair    of    atoms    changes    as    1 1 
distance    between    the    atoms    changes.       Initially,     at     large 
separations,      the     force     is     attractive,     becoming     stronger     as 
they    approach    one    another.    At    small    separations    the    force 
become     repulsive.  This        force     is        best     described     by     a 

potential     energy     function.  The     minimum     of     this     potential 

function  is  the  equilibrium  position  between  the  pair  of 
atoms.  Figure  10  shows  a  representative  plot  of  this 
potential    energy    for    a    pair    of    atoms. 


Repulsive 


Attractive 


Figure  10.   Representative  Potential  Energy  Plot. 

In  a  crystal  structure  this  potential  energy  function  is  a 
multi-body  problem.  The  result  of  these  extra  bodies  is 
that  at  a  greater  distance,  the  function  oscillates  about 
the  zero  energy  position  as  shown  in  Figure  11. 
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Figure  11.   Energy  Potential  of  a  !iulti-3ody  System. 


An  atom  has  a  kinetic  energy  that  is  characteristic  of 
the  thermal  energy  of  its  surroundings.  This  motion  is 
known  as  the  thermal  motion  of  the  atom. 

When  the  thermal  motion  and  the  interatomic  potential 
of  the  crystal  structure  are  combined,  the  atom's  notion 
becomes  a  vibration.  An  atom  oscillates  about  its  site, 
always  moving  near  its  lattice  position.  Add  a  diffusing 
atom,  in  an  interstitial  site  of  the  crystal,  with  its  own 
vibration  and  its  own  potential  energy  and  this  interstitial 
atom  will  gain  and  lose  energy  as  it  vibrates  between  the 
lattice  atoms.  Figure  12  shows  a  qualitative  graph  of 
energy  versus  time  of  an  atom's  history  while  in  an 
interstitial  site. 


Figure  12.   Time  History  of  an  Atom's  Energy. 

The  potential  energy  of  the  atoms  in  a  crystal  structure 
sets  up  an  energy  barrier  that  a  diffusing  atom  must 
overcome  before  it  can  move  to  one  of  its  neighboring  sites. 
This  energy  is  known  as  the  activation  energy  for  diffusion. 
Figure  13  shows  a  qualitative  graph  of  an  atom's  kinetic 
energy  with  respect  to  time,  with  the  activation  energy  for 
a  move  indicated. 


Time  of  jump  to  neighbor  site 

Figurel3.   History  of  Atom's  Energy  Noting  Atom's 
Activation  Energy. 
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D.   RANDOM  WALKS 

Early  in  the  study  of  diffusion  mechanisms,  t 1 
possibility  of  simulating  the  overall  motion  of  diffusion 
was  investigated.  The  method  devised  is  still  used  today 
and  is  known  as  a  Random  Walk  Process.  A  random  walk  in  one 
dimension  can  be  described  as  assigning  a  direction,  either 
forward  or  backward,  to  the  faces  of  a  coin.  When  the  coin 
is  flipped,  one  step* is  made  backward  or  forward  according 
to  the  face  of  the  coin  displayed.  Let  R  be  the  vector 
connecting  the  origin  to  the  final  position  of  the  atom. 
The  equation  of  R   is  then: 


Rn=r1+r2+...=    £ri 


(6) 


where  r  ^  arc  the  steps  representing  the  individual  jumps  of 
the  atom.  If  each  direction  of  motion  is  equally  likely  and 
individual  jumps  are  of  equal  magnitude,   then  the  magnitude 


of  Rn  is 


(Rn  )1/2=  ml/2 


(7) 


Now  let  a  number  of  atoms,  initially  at  x  =  0,  each  complete 
one  dimensional  random  walks.  Close  to  95%  of  the  atoms 
should    be    within    +    or    -    /i'       of     the    origin. 

S.       COMPUTER    SIMULATION    OF    DIFFUSION 

Shortly    after    the    development    of    the    electronic    computer 
Mayer     and    U 1 a  m     proposed     that     computers     be     used     to    simulate 
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physical  systems,  Metropolis  (1953), and  !'  i  n  r-  (1931) 
suggested  that  the  1 1  o  n  t  e  Carlo  method  of  simulation  ( t  h i  s 
net  hod  will  be  discussed  in  depth  below)  could  be  used  to 
simulate  the  random  walks  of  particles.  The  problem  w  a  s 
that,  at  the  time  of  King's  paper,  the  computational 
capabilities  of  the  machines  in  use  were  easily  exceeded  by 
Monte  Carlo  simulations.  By  1961,  computer  techno logy  was 
sufficiently  advanced  that  Flinn  and  HcManus  were  able  to 
complete  a  Monte  Carlo  simulation  of  solid  state  diffusion. 
However,  as  late  as  1970,  there  still  had  been  no  simulation 
of  the  diffusion  processes  based  on  the  random  walk  theory. 
By  1970,  it  was  evident  that  a  new  simulation  process  was 
required  since  the  analytical  processes  of  diffusion  were 
unable  to  solve  the  topical  diffusion  problems  evolving  in 
the  literature.  In  1971  Bennett  and  Alder  published  the 
results  of  their  simulation  of  diffusion  by  vacancies  and 
divacancy  random  walks.  Their  paper  marked  the  beginning  of 
the  random  walk  simulations  of  diffusion.  Through,  the 
1970's  and  198 0's,  the  field  of  Monte  Carlo  simulation  of 
the  diffusion  process  came  into  its  own.  Simulations  have 
been  used  to  test  the  correlation  factors  of  the  diffusion 
coefficient  (Guy  e  t  a 1 ,  1977;  Guy,  1973),  the  self  diffusion 
of  tracer  atoms  in  an  alloy  (Rakker,  1976),  and  many  other 
specific  diffusion  mechanisms  (.'lurch,   19  84). 
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F.'     CONCENTRATION    DEPENDENT    DIFFUSION    COEFFICIENT 

Concentration  dependent  diffusion  was'  apparent  early  in 
the  study  of  diffusion  and  the  diffusion  coefficient.  C . 
M  a  n  t  a  n  o  (19  3  3)  made  the  first  attempt  to  mathematically 
solve  for  the  experimental  results  that  showed  the 
concentration  dependence.  Using  the  experimental  results  of 
metal- metal  diffusion,  he  found  that  the  concentration 
dependent  diffusion  "coefficient  would  have  the  following 
form  : 


D(C)    =    (-l/2t)     (dx/dC)    /    xdC 


(8) 


This     equation     is     known     as      the     Matano     solution      for     the 
diffusion     coefficient. 

Bowker  and  King  (1978)  attempted  to  develop  a  Monte 
Carlo  simulation  to  determine  the  concentration  dependent 
diffusion  coefficient.  Their  first  attempt  used  a  two 
dimensional  lattice  gas  model  of  a  simple  square  plane  and 
their  results  appeared  to  correspond  with  the  experimental 
data  of  ox y gen-tun gston  diffusion.  However,  a  more 
detailed  investigation  of  the  simulation  led  to  discovery  of 
problems  in  the  data  collection  of  the  concentration  ( M u r c h , 
1981).  Hurch  (1980)  also  t  r  i e d  to  develop  a  simulation  that, 
could  provide  a  concentration  dependent  diffusion 
coefficient  (this  simulation  is  discussed  belo  w  )  ,  b  u  t  h  i  s 
simulation    showed    no    concentration    dependence. 
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A  recent  paper  by  A.J.  S  1  a  v  i  n  a  n  c!  ?  .  V  .  Under  hill 
(19  84)  attempted  to  prove  that  a  concentration  dependent 
diffusion  coefficient  could  not  be  found  u  s  i  n  ^  a  random  walk 
simulation  of  diffusion.  However,  a  paper  by  R.  Ghez  and 
U.S.  Lan^lois  (1966)  showed  that  a  concentration  dependent 
simulation  can  be  developed  if  care  is  taken  in  the 
collection  of  the  fluxes  at  the  midpoint  between  the  planes 
of     the    crystal. 


II.       OBJECTIVE 

The  primary  objective  of  this  thesis  is  to  investigate 
the  possibility  that  the  diffusion  coefficient  is  position 
dependent.  The  idea  that  such  a  dependence  night  exist  is 
not  new.  In  fact,  a  positional  dependence  lias  been  observed 
in  a  crystal  when  there  is  a  temperature  gradient  across  the 
structure.  It  has  always  been  accepted  that  w  he  n  t  hi  i  s 
gradient  is  absent,  and  the  temperature  is  constant 
throughout  the  structure,  the  diffusion  coefficient  would  be 
independent  of  the  position  of  the  diffusing  atom.  The 
idea  of  a  position  dependent  diffusion  coefficient  was 
proposed  by  Collins  (1985)  as  one  way  to  test  the  assumption 
that  the  partition  function,  in  Gibbsian  statistical 
mechanics,  does  not  divide  into  purely  positional  and 
velocity     parts. 

Though  proposed  early  in  the  history  of  simulation, 
diffusion  simulations  did  not  actually  begin  until  Bennett 
and  Alder  (1971)  developed  their  program.  They  were  studying 
persistence  effects  in  hard  sphere  gases  and  usee'  a  simple 
Monte  Carlo  simulation  to  find  the  vacancy  correlation 
factor  in  a  two-dimensional  fee  crystal.  The  first  actual 
simulation  of  three  dimensional  diffusion  was  done  by 
d  e  B  r  u  i  n  and  Hurch  (1973).  This  simulation  attempted  to 
develop     a     method     to     determine     the     diffusion     correlation 
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factor.        The    diffusion    correlation    factor    (f)    is    based 
the    Einstein    diffusion    coefficient    method    o  i."    determining    the 
coefficient.        Redefining    Equation     5 , 


D    =       n\Z    f/2t 


(9) 


the  diffusion  corellation  factor  can  be  defined  as 


f  =  <x2>/nx2 


(10) 


where  n  is  the  number  of  jumps  the  diffusing  particle  takes, 

2 
X     is  the  distance  of  each  jump,  and  <x>  is  the  mean- square 

displacement  of  the  diffusing  particle  after  time  t   (Le 

Claire,  1970).  Early  simulations  of  diffusion  dealt  with 

either  diffusion  of  a  vacancy,  or  self  diffusion,  and   all 

were  based  on   the  Einstein  diffusion  equation  (see  Equation 

5). 

The  first  reported  use  of  Pick's  First  Law  in  a 
diffusion  simulation  was  by  Guy  et  a  1 .  (19  77)  and  Guy  (1978) 
when  the  correlation  factor  was  calculated.  The  method  used 
in  these  investigations  was  to  place  an  initial 
concentration  of  vacancies  in  a  chemically  homogeneous 
square  planar  structure  and  then  allow  the  vacanc i e s  t o 
move  throughout  the  structure  b v  a  random  walk  process  for  a 
period    of    time.  This    method    was    found    to    have    problems     in 

that  they  had  assumed  that  the  atomic  flux  and  the  vacancy 
flux  across  the  plane  were  different,  and  set  up  the 
simulation     in     that     vein.         In     1979     Hurch    and     Thorn     a  n a  i  n 
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tried  to  simulate  the  tracer  correlation  factor  by  this 
method;  Their  method  was  to  set  up  a  three-dimension  a  1 
array  with  an  inexhaustible  reservoir  of  tracers  or 
vacancies  at  one  end  ( x  =  0 ) ,  and  at  the  other  end  (  x  =  1 ) 
remove  the  tracers  into  a  sink.  The  array  was  periodic  in 
the  y  and  z  directions.  From  this  simulation  one  could 
calculate  the  net  flux  across  a  plane,  and  from  this 
calculate  the  correlation  factor.  In  their  paper,  ^urch  and 
Thorn  reported  that  this  method  gave  as  accurate  answers  as 
the    Einstein    method    and    was    not    so    difficult    to    use. 

Murch  (19  SO)  also  used  the  Einstein  method  (Equation 
5)  to  test  the  effect  of  diffusing  atoms  concentration  on 
the  chemical  diffusion  coefficient.  This  simul.ation  was 
based  on  a  three-dimensional  simple  cubic  array.  Within  this 
array,  he  placed  two  reservoirs  of  atoms  eight  lattice  units 
in  from  the  x  boundaries,  each  with  a  separate  chemical 
potential.  He  then  allowed  the  atoms  to  diffuse  using  a 
nearest-neighbor  interaction,  interstitial  mechanism.  He 
concluded  from  these  investigations  that  the  concentration 
did    not    affect     the    diffusion    coefficient. 

Pick's  second  law  was  not  exploited  in  a  diffusion 
simulation  until  1978,  when  B  o  w  k  e  r  and  King  used  it  to 
study  the  diffusion  coefficent.  After  initializing  a  Lattice 
of  vacancies,  they  placed  a  concentration  of  defects  in  one 
half  of  the  lattice,  and  a  different  concentration  in  the 
other     half.        After     a     number     of     random     steps     of     these 


defects,  the  concentration  across  the  structure  was  counted 
plane  by  plane  and  a  concentration  profile  was  developed. 
From  this  profile  the  diffusion  coefficient  was  derived 
numerically.  Murch  (19  81)  used  the  second  la  w  n  e  t  h o  d  but 
with  the  external  reservoirs  that  he  and  Thorn  (1979) 
developed  for  their  first  law  investigation.  He  noted  that 
if  he  used  different  data  planes  to  develop  his 
concentration  profile'  he  found  different  results.  This  led 
him  to  believe  that  the  method  had  some  problems  with  the 
development  of  the  concentration  gradient  across  the 
structure. 

The  development  of  these  simulation  techniques  lias 
proven  that  the  diffusion  coefficient  can  be  studied  using  a 
digital  computer.  Each  of  these  simulations  used  the 
premise  that  the  coefficient  is  constant  throughout  the 
structure,  but  M urch  (1980)  suggested  that  these  techniques 
could  be  used  to  see  whether  the  diffusion  coefficient  has  a 
positional  dependence.  Another  limitation  to  all  of  these 
simulations  is  that  the  authors  used  either  the  Einstein 
method  to  determine  the  diffusion  coefficient  or  Pick's 
first  or  second  law.  The  simulations  were  never  analyzed 
with  more  than  one  technique. 

A  comparison  of  the  diffusion  coefficient  as  calculated 
by  Kick's  first  and  second  laws  will  be  the  primary  met  h  o  d 
to  determine  the  positional  dependence  in  this  thesis.  Tf  a 
positional  dependence  is  recorded,   a  method  other  than 
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Pick's  laws  will  be  necessary  to  determine  the  diffusion 
coefficient.  This  method  is  known  as  the  " F  o  k  k  e  r - P  1  a  n k ' 
equation    of    diffusion: 

3c/3t    =    -3     /3  x[A(x,  t)C(x,  t)  ]+    1/232    /3  x  2  [  D  (  x  ,  t  )C(  x  ,  t )  ]       (11) 

In  this  Equation  (11)  the  first  element  of  the  sun 
constitutes  a  drift  factor  where  A(x,t)  is  a  drift  velocity, 
and  the  second  element  is  the  diffusion  equation.  The 
equation  actually  is  giving  the  time  evolution  of  the 
probability  function,  C(x,t).  In  the  simulation  this 
probability  function  is  the  concentration.  If  the  time 
derivitive,  spatial  derivatives,  and  the  functional  form  of 
the  diffussion  coefficient  are  known,  the  diffusion 
coefficient    can    be    fit    to    the    simulation    data. 
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Ill .  COMPUTER  MODEL  AND  SIMULATION 

A.   PHYSICAL  MiODEL 

A  computer  simulation  model  must  mimic  a  physical  model 
or  experiment.  In  this  simulation  the  physical  model  is  a 
simple  face-centered-cubic  crystal  membrane  which  i  s 
infinite  in  the  x  and  y  axis  and  very  thin  (less  than  50 
lattice  planes)  along  the  z  axis.  The  mechanism  will  be 
interstitial  diffusion.  This  means  that  the  diffusing  atoms 
are  very  small  in  comparison  to  the  size  of  the  lattice 
atoms.  These  diffusing  atoms  start  in  an  inexhaustible 
reservoir  outside  the  input  plane  of  the  crystal, and  are 
completely  absent  from  the  volume  outside  the  other  "side  of 
the  crystal  at  time  equal  to  zero.  At  this  starting  time, 
the  atoms  begin  to  diffuse.  After  a  period  of  time  the 
distribution  of  diffusing  atoms  approach  a  steady  state;  in 
other  words  the  transient  effects  on  the  concentration 
profile  have  ceased.  At  this  time  the  flux  across  each  plane 
in  the  crystal  structure  is  checked  and  the  concentration 
gradient  is  recorded.  The  diffusion  coefficient  is  then 
found  from  these  values.  This  thesis  consists  of  a  computer 
simulation  of  this  physical  system,  which  w  a  s  designed  to 
study  the  diffusion  coefficient. 
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.       COMPUTER    SIMULATION 

The  physical  model  is  trans]  ate ci  in  such  a  way  that  t 
computer  can  follow  the  actio  n  s  t  hat  t  a k e  place.  ( I  )  i  .:  i  = 
known  as  computer  simulation  or  c  o  m  p  u  t  e  r  m  o  d  e  1  i  n ;  o  f  ; . 
physical  system.)  Computer  simulations  can  be  separat e  u 
into  two  families:  continuous-time  simulation  anr.  discrete- 
event  simulation.  Continuous-time  simulations  (also  known 
as  timestep  simulations)  simulate  events  w  h  i  c  h  c  a  n  1.  e 
described  by  simultaneous  equations.  The  equation s 
describing  events  are  allowed  to  proceed  for  a  pre- 
determinedperiod  of  time  (a  timestep),  and  then  information 
is  placed  back  into  the  equations  which  are  to  be  s o  1  v e  ■ ' 
numerically.  Discrete-evont      simulations      are      used      for 

physical  systems  which  cannot  be  described  by  a  mathematical 
function.  The  simulation  runs  from  event  to  event  rather 
then    advancing    a    number    of    actions    during    a    timestep. 

Simulations  can  also  be  divided  into  types  by  t  h  e 
methods  used  to  make  decisions.  If  the  simulation  u s c s 
r  a'ndom  numbers  to  make  decisions,  it  is  called  a  '. .  o  n  t  c  C  a  r  1  o 
simulation.  This  is  because  early  simulations  (those  prior 
to  computers)  used  devices,  sue::  as  dice,  to  choose  the 
random  numbers  that  were  used  to  make  the  decisions.  Once 
simulations  became  more  sophisticated,  -the  community  tried 
(in  vain)  to  shift  to  the  name  "stochastic  simulation',  but 
the  name  'Monte  Carlo'  ii a s  stayed  popular.  In  the  other  type 
of      simulation,       called      'deterministic      simulation',       all 
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decisions  are  made  by  rules  set  down  prior  to  the  actual 
running  of  the  simulation  program  (Harrison, 1985).  The 
simulation       used    in    this    thesis    is    of    the    Monte    Carlo    type. 

C.       RANDOM    WALK 

The  diffusion  events  are  simulated  by  a  method  called  a 
three  dimensional  random  walk  process.  This  process  is 
simple  in  form  and  can  easily  be  set  up.  Consider  a  t  w o 
dimensional  example:  Imagine  that  a  marble  is  on  a  Chinese 
Checkers  board.  The  marble  has  six  neighbor  sites  to  which 
it  can  be  moved.  A  random  walk  of  the  marbles  can  be  devised 
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here  each  direction  that  the  marble  can  move  is  assigned  a 


number  between  one  and  six.  A  die  is  rolled  and  the  marble 
is  moved  in  the  direction  indicated  by  the  number  on  the 
face  of  the  die.  This  is  an  example  of  an  uncorrelated  two- 
dimensional  random  walk.  Other  forms  of  random  walks  are 
the  % weighted'  and  ^correlated'  random  walk.  Returning  to 
the  Chinese  Checker  board  example,  if  you  move  two  marbles, 
one  that  is  black  and  one  that  is  white,  with  the  same  roll 
of  the  die,  that  is,  if  every  time  the  die  is  rolled  for  the 
white  marble  you  also  move  the  b  1  a  c  k  marble,  then  the  t  v;  o 
moves  are  connected  and  the  random  walk  is  correlated.  A  n 
e  x ample  of  a  weighted  random  walk  can  be  seen  by  t  a k i  n g  t  w o 
dice  and  assigning  one  number  to  each  of  the  four 
directions  you  want  to  have  a  lesser  chance  of  the  marble 
moving  to,  and  assigning  three  numbers  to  each  of  the  other 
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two  directions.  This  is  a  weighted  random  walk  in  the 
direction  of  interest.  The  random  walk  used  in  this  thesis 
is  a  Qon- weighted,  uncorrelated  random  walk. 

D.   COMPUTER  CODE 

The  code  for  the  computer  simulation  is  divided  into  two 
separate  programs,  x  D I F  S  E  T '  (for  diffusion  simulation  set-up 
program)  and  *  D I F  F  U  S  E  '  (for  diffusion  simulation  program). 
These  programs  combine  to  form  a  M onte  Carlo  simulation  of 
an  interstitial  mechanism  type  of  diffusion. 
1 .   Pi  f set 

The  program  'dif set'  can  be  separated  into  three 
parts,  each  of  which  intializes  a  table: 

TABLE  3 
TABLE  OF  DIFSET  ARRAYS 
Crystal  site  table 
Nearest-neighbor  table 
Site  vs  plane  table. 


a.   Crystal  Set  Up 

The  face-centered-cubic  crystal  used  in  this 
simulation  is  set  up  in  a  three  dimensional  array  called 
CRYST3(I,J,K).  The  t  h r  e e  dimensional  array  variables  must 
meet  certain  limitations.  These  limitations  are  that  the  I 
and  J  variables  have  to  be  even  numbers.  There  are  no 
limitations  on  the  K  variable.  The  reason  for  these 
limitations  is  that  the  1  and  J  axis  in  the  simulated  fee 


crystal  are  s  o  t  up  to  have  periodic  boundaries.  The 
initial  array  is  a  perfect  face-centered-cubic  crystal  with 
no  defects  or  vacancies.  An  array  position  designates  t  1 1  e 
center  of  a  site  and  gives  no  indication  of  the  size  of  an 
atom  that  might  occupy  the  site.  Atom  size  is  not  a  factor 
in  the  simulation  or  in  the  random  walk  used,  and  therefore 
is  not  included  in  the  set-up  beyond  the' assumption  that  the 
diffusing  atoms  are  "small  enough  (relative  to  the  lattice 
atoms  in  the  crystal)  to  allow  the  interstitial  mechanism  of 
diffusion  to  occur.  The  set-up  of  the  array  CRYST3(I,J,K) 
is  accomplished  by  placing  the  number  two  (2)  in  the  array's 
register  as  an  indication  that  a  lattice  atom  is  in  that 
array  position.  To  accomplish  this,  all  the  array  positions 
where  the  I  ,  J  a n-d  K  variables  add  to  an  odd  number  are 
assigned  the  number  two  (2).  Table  4  lists  examples  of  the 
initial  value  set  up  in  *  DIFSET"  for  the  lattice  positions. 


TABLE  4 


EXAMPLES  OF  ARRAY  ENTRIES  TO  THE  CRYST3  ARRAY 


CRYST3(3, 1,1) 
CRYST3(2,2, 1 ) 
CRYST3(2, 1,2) 
CRYST3(5,2,4) 


=  2 


? 


_  o 


The  interstitial  positions  of  the  crystal  are  set  up  in  the 
same  way.    Site  indices  that  add  to  an  even  number  are 
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assigned  a  zero  value  to  indicate  that  the  interstitial 
postitions  are  empty.  The  crystal's  first  plane 
interstitial  sites  are  filled  by  diffusing  atoms.  Therefore 
the  array  positions  (even,odd,l)  and  (odd,even,l)  are  given 
the  value  of  one  (1)  to  indicate  that  these  interstitial 
sites  are  filled  with  the  diffusing  atoms.  The  array  is  no  w 
set  and  an  example  plane  is  shown  in  Table  5. 

TABLE  5 

EXAMPLE  OF  THE  FIRST  TWO  PLANES  IN  A  4X4XK  CRYSTAL. 

CRYST3(I,J,K)    plane    #1 

2  12  1 
12  12 
2  12  1 
12    12 


CRYST3(I,J,K)    plane    #2 

0  2  0  2 

2  0  2  0 

0  2  0  2 

2  0  2  0 


The  three  dimensional  array  that  was  set  up  i  s 
difficult  for  the  computer  to  use,  so  a  one  dimensional 
array,  C  P.  Y  3  T  ( L )  ,  is  set  up  to  be  equivalent  to  the 
C R  Y S T 3 ( I , J , k  )  three  dimensional  array.  For  methodology  see 
code     i-n    Appendix    B , 
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b.   Nearest- Neighbor  Table 

The  next  major  step  of  the  *  D I F  S  H I '  progr  a  :.  i  s 
to  set  up  a  nearest-neighbor  table  that  will  be  usee,  to 
simplify  the  r  a  n  d o  m  walk  process.  To  review  the  diffusion 
by  interstitial  net  hod,  there  are  twelve  interstitial  sites 
which  are  nearest-neighbors  to  any  other  interstitial  site. 
Therefore,  a  diffusing  atom  must  move  into  one  of  these 
twelve  nearest-neighbor  sites.  The  array  C R R K B R  (  L  ,  R  )  (for 
crystal  nearest-neighbor  array)  is  set  up  such  that  the 
value  held  in  the  register  is  the  L  variable  of  another 
interstitial  site  in  the  CRYST(L)  array.  Thus  if  you  want 
to  know  the  twelve  nearest-neighbors  to  CRYST(13),  you 
would  look  the  it;  up  in  the  CRNNER(13,R)  where  the  twelve 
values  are  the  twelve  nearest-neighbor  site  L  variables. 
The  periodic  boundary  conditions  for  the  I  and  J  directions 
are  set  up  in  this  array.  The  boundaries  in  these 
directions  are  made  periodic  by  setting  up  the  neighbor 
t  a b 1 e  to  take  an  atom  that  leaves  t  h e  b  o  u  n  d  a  r y  of  t  h e 
crystal  on  the  left  side,  and  replace  it  on  the  right  side 
of  the  crystal  as  if  the  left  and  right  side  were  neighbor 
pianos.  The  edge  sites  of  the  K  axis  are  the  entry  an  d  e  :•;  i  t 
positions  of  diffusing  a  t  o  m s  and  as  such  have  special 
nearest-neighbor  table  values.  The  entry  edge  is  si  Tula  ted 
such  that  the  plane  is  in  contact  with  a  infinite  supply  of 
the  diffusing  atoms,  so  the  interstitial  sites  arc  always 
filled  with  the  diffusing  atoms.   Another  way  of  loo  icing  at 
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t  h  i  sa  s  s  u  m  p  t  i  o  n  physically  is:  if  a  diffusing  atom  leaves  a 
surface  interstitial  position,  another  diffusing  atom  takes 
the  interstitial  position  before  the  simulated  tin  estop  is 
completed.  Therefore,  the  neighbor  table  is  set:  up  to 
contain  only  the  four  neighbors  that  are  within  the  crystal. 
This  means  that  a  diffusing  atom  has  the  opportunity  to 
enter  the  crystal  only  one  third  of  t.he  time.  An  example  of 
this  table  is  located,  in  Appendix  F. 

c.   Site  vs  Plane  Table 

The  final  item  that  the  program  *  D I F  S  E  T ' 
develops  is  the  lattice  si  t  e  plane  table.  This  array 
ZPLN(L)  (representing  the  plane  §  in  --the  z  direction) 
correlates  the  site  number  L  to  the  plane  number  which  is 
contained  in  ZPLN(L). 

The  information  developed  in  '  D  I  F  S  E  T  '  is  then 
placed  into  an  output  file  to  be  used  as  input  to  the  other 
program  'DIFFUSE'. 

2 .   Diffuse 

%  D I F FU S E '  is  the  actual  simulation  program  that 
models  the  diffusion  process.  This  is  done  using  a  Monte 
Carlo  tinestep  simulation.  At  each  timestep,  the 
s i m u 1  a t i o n  looks  at  each  crystal  site.  If  it  contains  a 
diffusing  atom  the  simulation  then  looks  at  a  neighbor  site 
selected  at  random.  If  the  neighbor  site  is  empty  the 
diffusing  atom  is  moved  into  that  neighbor   site.   In  the 


35 


simulation,  the  process  is  done  with  two  tables  of  random 
numbers.  At  the  beginning  of  each  timestep,  a  random  number 
generator  is  called  to  make  a  table  of  random  numbers,  one 
for  each  site  in  the  C  R  Y  S  T  (  L  )  array.  The  random  number  is 
uniform,  1  to  L M AX,  where  L M A X  is  the  total  number  of  sites 
within  the  array.  The  sites  are  characterized  in  this  way, 
so  the  simulation  can  be  us'ed  later  to  study  diffusion  by 
the  vacancy  mechanism  with  little  modification.  Another 
table  of  random  numbers  uniform  1  to  12  is  set  up.  Again  the 
size  of  the  table  is  equal  to  the  number  of  sites  in  the 
array  CRYST(L).  As  a  time  saver,  during  the  first  few 
timesteps  the  number  of  random  numbers  generated  is  equal  to 
the  number  of  sites  in  one  plane  multiplied  by  the  timestep, 
and  only  those  first  planes  are  used  in  the  simulation. 

After  the  random  numbers  for  the  timestep  are 
assigned  to  their  tables,  the  simulation  calls  the  first 
random  number  (from  the  uniform  1  to  L M A X  set),  sets  that 
number  equal  to  L  and  looks  to  see  if  CRYST(L)  contains  a 
diffusing  atom.  If  the  site  does  not  contain  a  diffusing 
atom,  another  random  number  is  selected  from  the  first 
table.  Tf  there  is  a  diffusing  atom  in  the  site,  a  random 
number  is  selected  from  the  second  table.  This  number  Ls 
set  equal  to  R  and  the  nearest  neighbor  table  is  used  to 
find  the  value  in  CRNKBR(L,R).  This  value  is  then  sot  equal 
to  LI  and  if  CRYST(Ll)  equals  zero  the  diffusing  atom  is 
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moved  into  the  new  location.   This  process  continues   until 
all  L  M  A  X  numbers  have  been  tested. 

During  each  of  these  moves  a  number  of  factors 
required  for  the  numerical  solution  of  the  diffusion 
coefficient  are  calculated. 

a.   Concentration 

The  concentration  of  atoms  in  each  plane  is 
determined  at  the  beginning  of  each  timestep  by  looking  at 
every  site  in  the  array.  If  the  register  contains  a  one  (  a 
diffusing  atom),  then  the  concentration  in  that  plane  is 
increased  by  setting  K  =  ZPLN(L),  and  adding  one  to  the 
value  already  in  the  array  CONS(THSTP,K).  The  concentration 
is  then  divided  by  the  number  P L N M A X ,  which  is  the  total 
volume  of  a  single  plane.  This  puts  the  concentration  into 
units  of  particles  per  unit  volume.  These  units  are  needed 
for  use  in  Pick's  law.  This  value  is  stored  in  the  array 
CVOL(TMSTP,K) . 

b  .   F 1  u  x 

The  flux  density  through  each  plane  is 
determined  when  an  atom  is  being  moved  from  the  old  site  to 
its  new  site.  Since  the  flux  density  is  actually  between 
two  planes,  by  definition  the  area  between  plane  K ( 1 )  on  the 
left  and  K  (  2 )  on  t  li  e  right  will  be  considered  part  of  plane 
u  (  1  ) .  This  means  that  if  an  atom  diffuses  from  K  ( 1 )  into 
K ( 2 )     then     the     flux     of     K ( 1 )     has     1 / P  L  N  M A  X     added     to     it. 
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However,  if  the  atom  :n  o  v  e  s  froc  K  (  2  )  hack  to  [<  (  1  )  t  h  e  flu  :< 
of  ::  (  1  )  is  decreased  by  1  /  P  L  N  M  A  X  .  Thus  F  L  U  X  (  T  M  S  T  P  ,  K  ( -1 ) ) )  = 
FLUX(TKSTP,K(1))  +  1/PLMMAX  if  the  atom  moves  from  K(l)  to 
K(2)  and  FLUX(TMSTP,K(1))  =  FLUX(TMSTP,K(1))  -1/PLKiIAX  if  it 
moves   from  K  (  2  )  to  K  (  1 ) . 

c.  First  Derivative 

The  first  derivative  of  the  concentration  with 
respect  to  positioa  is  calculated  by  subtracting  the 
concentration  per  unit  volume  of  the  plane  after  the  plane 
where  the  derivative  is  required,  from  the  plane  preceding 
the  plane  where  the  derivative  is  required.  For  example 
DERC1(TMSTP,K)  =  CVOL(TMSTP  ,  !C  +  1  )  -  CVOL(TMSTP  ,K-  1 ) .  The 
units  for  the  concentration  gradient,  as  stated  above,  are 
particles  per  inter-planar  distance  to  the  fourth  power. 
With  these  units  the  c  h a  n  g  e  in  the  position  is  1  and  the 
mathematics  is  simplified. 

d .  Second  Derivative 

The  second  derivative  with  respect  to  the 
position  is  obtained  in  the  same  manner,  with  the  exception 
that  the  1st  derivitve  is  used  rather  than  the  planar 
concentration.  DE R C 2 ( T M S TP , K  )      =      DERC1(TMSTP,K  +  1)      - 

D  V,  RC1(TMSTP,K-1).  The  units  for  the  second  derivative  a  r  e 
particles    per     inter-planar    distance    to    the     fifth     power. 
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e .       Time    Derivative 

The  time  derivatives  are  derived  in  the  s ::  r.  e 
manner  as  the  positional,  derivative.  T  h  e  concentration  of 
the  plane  at  the  beginning  of  the  timestep  is  subtracted 
from  the  concentration  at  the  end;  DERT(T;IST?,K)  = 
CV0L(TMSTP+1,K)  -  CVOL(TMSTP,K).  The  units  for  this  factor 
are    particles    per    tin estep     per     inter-planar    distance    cube''. 

f  .       Averagin g 

After  a  pre-de ter m ined  number  of  timesteps  have 
been  completed,  the  arrays  containing  the  information  arc 
averaged.  This  gives  another  set  of  arrays  that  are  used  to 
derive  'the  diffusion  coefficient.  The  single  timestep 
arrays  are  then  reset  to  zero,  and  another  set  of  timesteps 
arebegun. 

g.       Diffusion    Coefficient 

In  this  simulation  the  diffusion  coefficient,  is 
derive  d  b  y  two  methods  so  that  t  h  e  results  can  be  cor:  p  a  r  e  d 
and  the  positional  dependence  of  the  coefficient  may  be  more 
easily  seen.  The  first  method  uses  Tick's  first  law. 
Having  found  the  average  flux  over  a  number  of  timesteps, 
AFLUX  (  ATi  ISTP  ,  K  )  ,  and  the  average  first  derivative, 
ADERC1(ATMSTP,K) ,  the  diffusion  coefficient  can  be  found  by 
dividing  the  average  flux  by  the  negative  of  the  derivative, 
DIFC01(ATMSTP,K)=  - AFLUX( ATMSTP , K) / ADERC 1 ( ATMSTP  ,  K)  .  The 
diffusion    coefficient     is     not     determined     for     every     timestep 
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because,  at  times,  t.h  e  first  derivative  for  a  sin'  5 
tin es ten  is  zero  and  the  resulting  diffusion  coefficient  is 
infinite.  This  result  w  o  u  1  d  skew  the  averaged  ri  i  f  f  u  s  i  o  n 
coefficient  . 

The  second  method  of  determining  the  diffusion 
coefficient  is  based  on  Pick's  second  law  of  diffusion. 
More  specifically,  the      second      law      holds      t  h  a  t        the 

diffusion  coefficient-  is  constant  with  respect  to  position, 
9  c  /  3  t  =  D  8ic/  9x*~.  The  diffusion  coefficient  is  then 
determined  by  dividing  the  time  derivative  of  the 
concentration  by  the  second  positional  derivative, 
DIFC02(ATMSTP,K)  =  AjQL&BiT  (  A  T  M  S  T  P  ,  K  )  /  A  D  E  R  C  2  (  A  T  M  S  T  P  ,  K  )  )  . 
A v  e  r  a  g  e  s  are  used  here  for  t  h  e  same  reason  as  described  for 
the     first     method. 


hi) 


tv.   nr.-VuLTS 

A.       PROOF    OF    SIMULATION'S    FFYSICAL    EQUIVALENCE 

The  proof  of  the  simulation's  ability  to  no  del  the 
physical  world  can  be  seen  in  three  ways.  First  the 
simulation  adheres  to  the  conservation  of  matter.  Second, 
the  diffusion  coefficient  that  is  developed  from  the 
simulation's  data  lies  within  the  limits  of  the 
experimentally  found  values  for  ■  the  diffusion  coefficient. 
Third, the  time  average  of  the  concentration  profile 
approaches  the  infinite  time  result  predicted  for  the 
concentration  . 

1 .  Conservation    o f    matter 

The  conservation  of  matter  and  Fick's  second  lav; 
(Equation  3)  are  equivalent.  Therefore,  if  Equation  3  holds 
in  the  simulation  then  the  conservation  of  matter  holds  by 
definition.  To  test  that  this  requirement  is  mot  at  all 
times  during  the  simulation,  the  first  derivative  of  the 
flux  density  was  determined  and  compared  to  the  time 
derivative  of  the  concentration  during  one  run  of  the 
simulation.  It  was  found  that  the  conservation  of  matter 
held    for    all    time    in    the    simulation. 

2 .  Comparison    o  f    t  he    Diffusion    Coefficient 

The  diffusion  coefficient  developed  from  the 
simulation    data    is    in    units    of    inter-planar    distance    squared 
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per  tineste?.  T  n  o  r  d  o  r  L  o  compare  fh  e  diffusion 
coi'.'fic  iont  round  "by  experiment  w  i  t  h  the  coefficient 
developed  from  the  simulation,  the  lattice  so  a c  i  n  g  of  I 
base  material,  the  tcr'Derature  of  the  system,  the  mas  s  o  f 
the  diffusing  particle,  and  the  force  constant  of  the 
harmonic  oscillator  approximation  to  the  potential  function 
between  the  diffusing  atom  and  the  1)  a  s  e  material  m  u  s  t  b  e 
provided.  As  an  ex  a*m  pie,  the  diffusion  coefficient  for 
carbon  diffusing  through  fee  iron  at  910  degrees  Celsius  has 
been  experimentally  determined  to  be  1.91  x  1 0 ~ '  c  m  / s 
(A  sic  eland,  19  84).   The  lattice  unit  for  fee  iron  is  3.589  x 

_  o 

10  cm  (Ask eland,  1 9  S  4 )  ,  the  force  constant  for  carbon- 
iron  is  in  the  range  2.0 5  x  10°  to  17.73  x  10'  dynes/cm 
(A.J.  Gordon,  1972),  and  the  atomic  w eight  of  carbon  is 
12.011  gr/mol.  From  these  constants,  the  experimentally 
determined      diffusion      coefficient  listed      above       is 

transformed  into  a  r  a  n  g  e  of  values  that  fall  between  0.23 
and    0.61     inter-planar    distance    squared    per    timestep.  '    - 

diffusion  coefficient,  values  determined  by  the  t  i  i 
averaging  over  10C1  timesteps,  fall  between  0.19  and  0 .  3 
inter- planar  distance  squared  per  timestep.  The  values  for 
the  coefficient  determined  b  y  t  h  e  s  i  m ulation's  t  i m  e 
a  veraging  over  1  0  C  0  timesteps  is  b  e  t  v:  e  e  n  0.2  6  a  n  d  0  .  4  2 
i  n  t  e  r - p 1 a n  a  r  distance  squared  per  timestep.  T h is  in  d i  c  a  t  e  s 
that  the  simulation  is  determining  the  diffusion  coefficient 
w i  t  h  i  n    acceptable    limits. 
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3.   Concentration  Profile' 


n  e  concentration  orofile  a  c  r  o  s 
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proceedtoa  1  i  n  e  a  r  slope  across  the  planes  of  t  h e  c  r  y  s  t a  I 
structure.  The  simulation  *  DIFFUSE'  begins  to  approach  the 
correct  forn  at  the  1-1000  tine  average  of  the  concentration 
(see  figure  14).  It  then  proceeds  to  oscillate  about  the 
infinite  time  profile  result.  (see  figure  15-17).  This 
oscillation  about  the  infinite  tine  line  indicates  that  the 
simulation  reaches  an  approximate  steady -state 
configuration.  To  test  this  'hypothesis,  the  19001-20000 
time  average  profile  was  statistically  checked  against  the 
linear  infinite  fit.  The  statistical  method  used  was  a 
linear  regression  of  the  data  to  the  infinite  result.  The 
regression  results  gave  a  probability  of  0 . 9 9 9  that  the  data 
does  fit  the  linear  slope  of  the  infinite  tine  profile.  The 
sane  approximate  result  was  found  for  the  14901-1500 0  a  n  d 
the  60901-61000  tine  average  concentration  profiles.  T  h c 
sa n e regression  was  fit,  using  the  data  averaged  over  19901- 
20  0  00timesteps,  f  o  u  n d  that  the  probability  t  h at  it  s 
concentration  profile  could  be  fit  to  the  infinite  time 
profile     w as     0.953. 

".j  oca  use  of  the  concentration  profile's  closeness  of 
fit  to  the  infinite  time  profile,  the  apparent  derivation 
of  the  actual  diffusion  coefficient  by  the  simulation,  and 
the    simulation's    ability    to    obey    the    conservation    of    matter, 


f,   o 

<-*  -J 


CONCENTRATION  PROFILE 


—4  ~ 

05 

\      '• 

Vy  ! 
\ V 

'.       \  \ 
;         \ 
V 

\  \    ' 

LEGEND 
INFINITE  TIME  PROFILE 

-• 

o~ 

•   TIMESTEP  1000  PROFILE 

CO 

VOL) 

.7          0 

1 

\     \  ; 
\     \ 

\  ;\ 
\  •  \ 

UNIT 

.6          0 

•   \      > 

;    \ 

;        \ 

s      \ 

PART/ 

.5          0 

\ 
\ 

N 

\ 
\ 

-    : 

0          V 

)  NOIJ 

5 

,       \; 
\        \ 

\     :\ 
N    •  \ 
\ ;    \ 

1 

s 

\     \ 

\       1   \ 

:once 

.2          0 

\          \ 
"  s         \ 

!             N 

\ 

\ 

N 

\  \        • 

\  \ : 

o 

o 
d 

1 

~" r* 

i    ■ 

i 

1   ^^^ 

V 

• 

■     "T 

12      3      4      5      6      7 


8      9     10     11     12    13    14    15    16 
PLANES 
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it  is  reasonable  to  assert  that  the  simulation  is  acting  as 
as  an  approximate  physical  model  of  the  interstitial 
diffusion  mechanism  through  a  fee  crystal. 

: .   CALCULATION  OF  THE  DIFFUSION  COEFFICIENT 

1  .   Tick's  First  L a v;  Method 

% DIFFUSE'   initially   develops   the   diffusion 


coefficient  by  using  the  F  i  c  k ' s  first  law  method. 
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m  e  t  h  o  d  takes  the  time  a  v  e  r  a  °  e  flux  and  divides  it  b  v  the 
time  average  first  derivative  of  the  concentration.  3 y  this 
method,  the  diffusion  coefficient  shows  what  appears  to  be- 
an oscillation  across  the  planes  (see  Figures  18  c.  19).  The 
cause  of  these  oscillations  could  be  the  result  of  either 
noise  in  the  simulation,  the  data  fitting  a  horizontal  line 
across  the  plane,  or-  a  positional  dependence  on  the 
diffusion  coefficient.  To  test  whether  the  data  fit  a 
'Horizontal  line  across  the  plane,  a  linear  regression  was 
fit  to  the  data.  The  results  of  the  fit  are  plotted  in 
Figures  13  3  19.  The  linear  correlation  factor  for  the  time 
average  gives  a  probabillity  of  20  -  30  per  cent  that  the 
data  fits  a  line  of  constant  slope  and  less  then  3  0 
percent  probability  that  the  data  fits  a  h orizontal  line 
across  the  planes. 

To  investigate  the  chance  that  the  oscillation  was 
noise  and  not  a  positional  dependence  on  the  diffusion 
coefficient,   the  diffusion  coeficient  data  as  a  function  of 
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Figure    18 
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It  is  a  p  p  a  r  e  n  t  f  r  o  ;r  the  t  v  o  t  a  b  1  e  s  t  li  at  t  h  e  results  o  f  t  '  i  e 
t  w  o  methods  are  drastically  different,  i'nli'c  t  i  i  e  \:  i  r 
r.\  e  t  h  o  d  data  the  second  la  w  data  is  not  equally  d  i  s  t  r  i  b  u  t  e  c 
a  b  o  u  ti  t  s  nean  a  n  d  a  1  1  t  h  e  data  are  n  o  t  w  i  t  h  i  n  t  w  o  standard 
deviations  of  the  rr.  e  a  n .  The  means  of  the  two  m  e  t  h  o  d  s  a  r  e 
also  net  w i  t  h  i  n  two  standard  deviations  of  each  o  t  h e  r . 
Therefore,  the  diffusion  coefficient  nay  be  dependent  on  the 
position  of  the  diffusing  atom.  The  next  step  in  the  test 
of  the  hypothesis  is  to  find  the  first  derivative  of  the 
diffusion  coefficient  and  compare  the  time  derivative  of  the 
concentration  and  the  result  of  Pick's  second  lav  with  the 
n  on- linear  diffusion  coefficient  (Equation  3a).  The  results 
of  this  comparison  can  be  seen  in  Tables  8  and  9. 

TABLE  OF  SHORT  TIME-AVERAGED  TIME  DERIVATIVES 


Plane 


Lime  Derivative 

.  0  0  0  5  6 

-.00056 

-.00222 

-.(.'01  11 

.00000 


F i c k ' s  Seco n d  L a w 


-.02396 
-.00320 

-.ooi9(; 

-.00481 

-.00  173 


iOTE:  Table  of  time  derivatives  from  simulation  and  derived 
by  Equation  3  a  us  in  3  the  averages  over  the  19  901-2 
tircesteps. 
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TABLE  OF  LOPG  TIKE-AVERAGED  TIIIE 


Plane 
3 
4 
5 
6 
7 


T  L m e  ;jer  L v a t i v e 
.00017 

.00000 
-.00006 
-.00011 
-.00  0  1 1 


.I\  ATI'*  -: 


i  c k '  s  Second  La' 


i   •-  ">  a 

W     \J    *m   ~J   -J     , 

-.('0142 

-.000  2  7 

.00104 

-.00107 


MOTE:  Table  of  time  derivatives  from  simulation  end  derive: 
by  Equation  3a  using  the  averages  over  the  1 9  001—2' 
timesteps . 

The  results  of  this  con  p  arisen  show  that  the  non-  linear  for:1, 
of  Pick's  second  lav.'  does  not  agree  wit  h  t  he  si  m  u  1  a  t  i  o  n 
results  of  the  tine  derivative  of  the  c  o  n  c  e n  t  r  a t  i  o  n .  This 
suggests  that  Pick's  first  law,  taken  literally,  may  not 
be  the  best  method  to  fit  the  simulation  data  to  a  diffusion, 
coefficient.  It  nay  be  the  diffusion  coefficient  should  be 
within  the  partial  derivative. 
3.   Fokker -Plank  he t hod 


The  last  method  used  to  search  for  a  positional 
dependence  was  to  attempt  to  fit  a  function  to  the  diffusion 
coefficient  and  solve  t  h  e  F o  k  k  e  r - P 1  a  n k  e  q  u a  t  i  o  n  ( E  q  u a t i o  n 
13)  with  this  function  and  the  results  of  the  si n u 1  a  t i o n . 
T  h  e  folio  w  i  .  n  g  functional  for  m  w  as  assumed,  Id  a  s  e  d  o  n  t  h  e 
diffusion  coefficient  resulting  from  the  Pick's  first  la  w 


me t  hod 


D(xft)  =  A  +  iisin(Cx  +  Dt  +  E) 


(  1  2 ) 


Using  this  equation  and  the  results  fro™  the  simulation,  the 
Fokker-Plank  equation  takes  the  following  form: 

3C/8t  =  Bsin(Cx+Dt+E)[0.5C2(CV0L)+C(DERCl)=.5(DERC2)  ]  +  A/2 

(13) 

The  IliSL  library  subroutine  ZXSSO  was  used  to  evaluate  the 
Fokker-Plank  equation.  Z  X  S  S  Q  finds  the  miniraun  of  the  s  u  r. 
of  squares  of  M  functions  in  M  variables  using  a  finite 
difference  Levenberg-Marquardt  algorithm.  In  the  simulation 
the  N  variables  are  the  A ,  B ,  C ,  D ,  and  E  variables  of  the 
diffusion  coefficient  function,  and  the  i I  functions  are 
Equation  13  evaluated  at  different  tin  estops  from  the 
simulation. 

It  was  determined  that  the  best  fit  for  the 
coefficients  A ,  B ,  C ,  D  and  E  was  w h en  they  all  e  q  u a  1  e  d 
zero.  This  result  indicated  that  the  oscillations  in  the 
diffusion  coefficient  './ere  probably  based  on  the  numerical 
noise  in  the  system. 


5  5 


V  .   CC  ..LL.l.  o  iC1..  :-i  A  .  .j  ..i-.^0...  ■ ilw,.- 

The  data  fro-  the  simulation  gave  no  c  1  ear 
determination,  one  way  or  the  other,  about  the  positional 
dependence  of  the  diffusion  coefficient.  The  simulation 
should  be  further  improved  to  nake  a  specific  determination 
of  the  positional  dependence  of  the  diffusion  coefficient. 
The  following  improvements  to  the  simulation  should  give 
it  the  ability  to  make  the  determination  of  the  dependence. 

First,  the  time  averaging  of  the  flux,  concentration, 
and  derivatives  may  have  averaged  out  any  real  positional 
dependence.  Therefore,  the  number  of  tiiuesteps  the 
simulation  averages  over  should  be  shortened,  or  deleted 
altogether.  Instead,  the  few  times  that  the  derivatives  are 
zero  should  be  ignored  and  the  diffusion  coefficient  should 
be  inspected  at  all  other  times  and  planes. 

Second,  the  size  of  the  array  should  be  increased.  The 
size  used  in  this  thesis,  6  X  6  X  1 6 ,  m ay  not  h  a v e  b e  o  n 
large  enough.  It  appears  that  with  this  size  a  single 
p  a  r  t  i.  c  1  e  change  in  a  plane  leads  to  large  changes  in  the 
parameters  determined  by  the  simulation.  The  array  should  be 


increaseu  in  size 


:    to    at    least    a     14    X    14    X    5  0    a  r  r  a  y • 


h  l 


size  should  make  a  single  particle  c i  i  a  n  g  e  in  a  plane  equal  a 
one  percent  change  in  the  determined  factors  of  the 
simulation  . 


r>5 


Third,  the      functional      form      of      the      d  i  f  f  u  s  i  a  n 

coefficient,  and  the  initial  values  for  t  h e  n n a  1 y s i  s  o C  t '  e 
u  n  h.  o  w  n  variables  cor  the  F  o  k  k  e  r  -  ?  1  a  n  k  equation  no  e  c'i  t  o  b  e 
investigated     further. 

A  major  result  of  this  thesis  is  its  improvement  in  the 
boo keeping  ability  of  diffusion  simulations.  *  D I F FUS E '  uses 
arrays  to  keep  track  of  values  that  are  needed  to  simulate  n 
more  sophistacated  m o d e 1  of  diffusion.  These  a r r s y s  k e  | 
track  of  the  nearest-neighbors  for  the  diffusing  a  tors,  and 
the  planes  that  the  atons  are  in.  This  a  1  1  o  i;  s  the 
simulation  to  be  written  more  simply  than  earlier 
simulations  and  allows  the  more  sophisticated  nodclin  '/•  o  f 
the     face-centered-cubic    crystals    to    be    sir;,  ulated. 

The    simulation's    improvements    to    the    method    by    which    the 
diffusion    process    is    modeled,     along    with    these    suggested 
improvements    to    the    simulation,     should    allow    the    simulation 
to     determine     if     the     d i  f  f  u  s  i o  n     coefficient     has    a     positional 
dependence . 
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VARIABLE  GLOSSARY 


The  following  glossary  of  variable  names  is  set  up  to 
aid  in  u  s  i  n  3  the  computer  programs  *  D I F S E T '  and  *  D I F F U S "  ' . 
The  glossary  consists  of  three  columns,  the  name  of  the 
variable,  the  type  of  variable  and  t  h e  p h y  s  i  c  a  1  m c  a  n  i  n  3  of 
the  variable. 


VARIABLE    NMVM1 


VARIABLE    TYPE 


:•>  • "  V  c  T  C   '•  T   ' :  7  '   ' '  T  ' "  i 


ACONS 


2-B    ARRAY 


coeceetrati 

EACH    PEA!"  I!  . 


AC  VOL 


2-B    ARRAY 


I  i  -    ■  r,    •   rr; 
. ■.  *   .-    .  n  u  i. 

CONCENTRATIOI 

PER   UNIT   VOL.    AT 

r  \  rT]     dt   <  •■  - 

j.j  r.l,  ;1       t    UrH<  u  < 


AD  FAR  CI 


2-D    ARRAY 


AVERAGE  1ST 
D  E  R  I  V  A  Y  I  V  E  I  A 
TEE  CONCRETE A- 
T  I  0  A       A  T       A  A  C 

RLAAE. 


ADERC2 


2-D    ARRAY 


AVERAGE  2  A  A 
DERIVATIVE  AT 
EACH    PEAAE. 


ADERT 


2-D    ARRAY 


AVERAGE       TIME 
DERIVATI  V E     OF 

1  :.  1/       U  Vj  :•!  L  l.  i>  i  n  ft  - 

TIC  A       AT       RACE 

PLAEE. 


A  FLUX 


2-D    ARRAY 


AVERAGE    EL A  A    OF 
EA'CA    PEAAE. 


ATMSTP 


ARGUMENT 


USED       AS       1ST 
ARGUMENT    IE    AAA 

TEE      A  V E  R  A  C  A . 


ARAAYS. 


ATI 


PARAMETI 


r  t  r:  - 


TI A  ESTEPS  A 
TEE  A  V  A  A  A  A 
ARRAYS. 


n 


AVGT 


PARAMETER 


t«   U   i!  B  r.   R 

TTMFCTFPC      v  r 

ft    V     [ju  ft  I.T   lj  ij  vj     V      ,  .    [  .      . 


a  \;  c  T  \f  1 

ft    V    \j    {    .  .    1 


PARAMETER 


PARAMETER     A 
TO    A OLA    EAAA 
T  T  A  A  S  T  !■  A  A 
MEMORY. 


]  i  i 
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VARIABLE  '.' r. 


V '  :IA?,LE  TYPE 


physical  feafin 


2-D    ARRAY 


:■  ncentra'  I   - 
for   eacf   pla 

Alio    TIliESTEP. 


CVOL 


2-D    ARRAY 


concentration 
per  unit  vol  at 
each   pl  art i 
tiiiestep. 


1 1<  i 


r  n  *  —  7  n 

k.  i\ . . .   .  ■  ■ 


O      p       a  n  n   i  V 
Z  —  !J       A  K  n  :l  I 


» ■  T7  inrrT      JTTr;!TiriTi 

i,  ;j.'l.....;  1  -:.  ,.  1':.-   >  l\ 

TABLE. 


CRY ST 3 


\  n  ~  \  v 


THREE  DI.M  :  ;£  10"- 

AL  ARRAY  OF  FCC 
CRYSTAL. 


p  n  y  C  T 


1-D  ARRAY 


oke-difefsional 
equivale::t   tc 

CRYST3. 


n  r  P  ^  i 


O  P  A    T)  p    A    V 

Z.  —  ly      rt.Jwil 


1ST  derivative 

OF    THE    COFCEF- 

T1  r    \  m  T  A  *  1         t  i — '       *~~    \   f^  T "" 

i  ;   ail  U  .-.     i\  i      Hi  u  V*  1 

PLATE. 


DERC2 


o     r,     a  n  p  a  y 


2TD     DERIVATIVE 
AT    EACH    PLAKE. 


nrDT 


2-D    ARRAY 


TIFF    DEIVATIVE 
AT    FACF    FLAFF. 


DIFC01 


-D    ARRAY 


DIFFUSICi:    COEF- 
FICIENT   AT    FACT 
FLAFF ,        DETER- 

-  j  ->  c  7     I  ii 


DIFC0  2 


_  P       A  T1  p   A 

.  —  J         iilw.il 


diffusion  coef- 
ficient deter- 
mined fy  fick's 
second  la::. 


FLUX 


2-0    AR  F\Y 


FFFF    ACROSS    FACT 


D  T     A  V  TT        A 

1         '_!    .  1 


TT 


.1-  . 


6  o 


variarl 


VA .  I  V-BLE   TV 


-v 


Y  S I C  A  L 


[    • 


i::a::i 


i . . ._  i  .. 


ix 

~  t  n  ■ 


THAN    IMAX. 


ISEED1 


raxdo 

c  —  — -n 
j  i  j  i  j  J  . 


iaa^aa 


DAP  ^HITTrn 


:.a:;  do: 

SEED. 


JMAX 


PARAMETER 


\'TT'  ;  p  -  -     ni7  •  13  I    \  vr> 

T  "'     V     1  T  r  "~  rr:  '  {    7 


JAAX1 


PARAMETER 


HA  A      AAA  A 
JAAA. 


aiibL.i  1.. .  i 


ARGUMENT    TO    / 

2-DIAAASIOAAL 


A  R  R  A  Y  S 


p 
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5b.ii5 


;-■'■,  v 

1V.J  lll.i 


NUMBER    OA 
I A    TAE    Z 


i  T     AM"; 


n 


■i  LA..  . 


K  1  .  j\  A  I 


OAE      LASS 


aaa: 


ARGUMENT 


A  A  G  U 


rt  i 


1  b 

CAYST. 


PT 


AAA  A 


p  ,    „   A  >  .♦  IJ  rp  j;  -p 


A  A  A  A  A  A    0  A    S  I  T  "  A 
IN   CRYST. 


LMAX1 


p   IT)    \     ■ TT  V  r> 

*    ai;a.._:iji> 


LMAX-PLXMAX. 

gaa    aaa:" a   a::: 

AAA  A       ALL       T- 
SITES     IX    CE'i 


lrd; 


,  -,  r  ...... ..  . 

.'...>  J  L  :   ...  .    i 


AAA  A. 


;r 


1      A 


^ 


AAAAYS 


VAR] '     L£ 


[apfle 


DTT 


.    1.  I    .  I  .  .     . 


'-1 


rL  .  ill/ 1.  1 


L    i:i  I   J  .... 

■  z:;t  n  . 
to  n:: 
neighbor. 


"•■'I  - 


.. . 


:iavgt 


PARAMETER 


largest    ;;u:  . 

OF    TllV.    A  VET.  AC 
TIMESTEPS. 


:sti 


D 


LARGEST     NUM! 

OF    TIMESTEP:    . 


OUTPUT 


n  tv-  ;  i 


LOGICAL 


1-D    ARRAY 


DETERMINES 
AMOUNT    OF   OUTP1  ^ 


r^     TT     -i                   T)     j-i      ,'~.     f 

ill'.            J.      1 ..     ..'   '  - 

GIVES. 

ARRAY    OF    FA- 

r. rs  *  ' 

FF. .FURS    USED 

FOR 

PICKING      C". 

. 

SITE        TO 

"■" 

LOOKED    AT. 

RDM2 


1-D    ARRAY 


ARRAY    OF    RANDOM 
NUMBERS    USFF    I  I 

D  T  r  V         V  T?   a   p   r  C  T 

NEIGHBOR    TO    FF 
LOOKED    AT. 


TFSTP 


ARGUMENT 


A  R  G  U 

T  H  F 


r  '•■  T 


zi 


A      i      _; 


1  /-  ^ 


J  i  'j    .  v 


AFF    CERT. 


i    n  ~  -  r  2 


Tl 


i  A  n  a  •  t  F  T  F  T? 


FIRST     TIMESTLP 
F  U  U  F  E  F  OF 

PROGRAM. 


ZPLN 


1-D    ARRAY 


::  i     it  ft  1         > .   l         .. •   I    i   . .  ;  i 

I'C      n  t     »  :  i  I •  c 

>  ."'      ,.    ij  .  i  . ,  i    j  • 
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APFEKDIX 


\    -     T  "  c    r  r-.    » 


IFSET'    CODE    LISTING 
The      following      is      the     computer     c  o  ci  c      for      the      set-ur 
program     *  D  I  F  S  E  T  ' .        This     code     develops     the     input     for     the 
program    ^DIFFUSE1 . 


r   n 
')  .J 


C  'DIFSET' 

C  MARK  R.  POLNASZEK 

C  JUNE  1986 

C 

c 

C    'DIFSET'  WAS  DESIGNED  TO  INTIALIZE  THREE  ARRAYS  THAT  ARE  NEEDED 

C    IN  THE  PROGRAM  'DIFFUSE'.   THESE  ARRAYS  CRYST(L),  CRNNBR(L,R), 

C    AND  ZPLN(L)  ARE  USED  IN  'DIFFUSE'  TO  SPEED  UP  THE  CALCULATIONS 

C    OF  THE  DIFFUSION  PROCESS. 

C 

C 

C    THE  FOLLOWING  DEFINES  THE  VARIABLES  OF  THE  PROGRAM  AS  REAL  OR 

C    INTEGER  VALUES,  AND  DIMENSIONS  THE  ARRAYS'  MEMORY  SIZE. 

C 

C 

00020   INTEGER  PLNMAX,CRYST(  2000 ) ,  CRYST3J 10 ,10,20 ) ,CRNNBR( 2000 ,12 ) ,R, 

1  ZPLN(2000),  ISEED1,  ISEED2,  FLUXI2000),  CONS( 2000 ) ,MTMSTP, 

2  T1,AVGTM,PLSTM, OUTPUT 

C 

c 

C   THE  FOLLOWING  READS  IN  THE  INPUT  TO  THE  PROGRAM  AND  ECHOS  IT 

C    BACK  TO  AN  OUTPUT  FILE. 

C 

C 

READ  (5,  *)IMAX,  JMAX,  KMAX,MTMSTP,T1 ,ISEED1,ISEED2,AVGTM,PLSTM, 
1  OUTPUT 

WRITE  (6,26)  IMAX,  JMAX,  KMAX,MTMSTP,T1 ,AVGTM 
00026   FORMAT  1 'l',6I8) 

WRITE  (6,28)  ISEED1,ISEED2 
00028   FORMAT  ( '1', 2116) 
C 
C 

C    THE  CODE  BELOW  INITIALIZES  CONSTANTS  THAT  ARE  USED  THROUGH  OUT 
C    THE  PROGRAM. 
C 
C 

00040   PLNMAX  =  IMAX  *  JMAX 
00050   LMAX  =  KMAX  *  PLNMAX 

PRINT  52,  PLNMAX,  LMAX 
00052   FORMAT  ( 'l',2I8) 
C 
C 

C    THE  FOLOWING  DO  LOOPS  SET  UP  THE  FIRST  ARRAY  ZPLN(L).   THIS  ARRAY 
C    HOLDS  THE  VALUE  OF  THE  PLANE  NUMBER  FOR  THE  GIVEN  SITE  NUMBER, L. 
C 
C 

DO  3192  L  =  1,LMAX 
ZPLN(L)  =  (INT  (  L/PLNMAX  )  HI 
DO  3191  N  =  PLNMAX, LMAX, PLNMAX 
ZPLN(NI=  N/PLNMAX 

03191  CONTINUE 

03192  CONTINUE 
C 

C 

C    THE  FOLOWING  DO  LOOPS  ARE  USED  TO  SET  UP  THE  THREE  DIMENSIONAL 

C    ARRAY  CRYST3(I,J,K).   THIS  ARRAY  HOLDS  THE  VALUE  OF  THE  TYPE 

C    OF  ATOM  IN  THE  CRYSTAL  SITE  OR  THE  FACT  THE  SITE  IS  EMPTY. 

C 

C 

00060   DO  102  K  =  2, KMAX, 2 

00062      DO  70  J  =  1, JMAX, 2 

00064         DO  68  I  =  1, IMAX, 2 

00066  CRYST3(I,J,K)  =  0 

00068         CONTINUE 

00070      CONTINUE 

00072      DO  80  J  =  2, JMAX, 2 

00074         DO  78  I  =  2, IMAX, 2 
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00076 
00078 
00080 
00082 
00084 
00086 
00088 
00090 
00092 
00094 
00096 
00098 
00100 
00102 
00104 
00106 
00108 
00110 
00112 
00114 
00116 
00118 
00120 
00122 
00123 
00124 
00125 
00126 
00128 
00130 
00132 
00134 
00136 
00138 
00140 
00142 
00144 
00146 
00148 
00150 
00152 
00154 
00156 
00158 
00160 
00162 
00164 
00170 
00172 
C 

c 
c 
c 
c 
c 
c 
c 

00174 

00176 

00178 

00180 

00182 

00184 

00186 

00188 

00210 

C 

C 


CRYST3(I,J,K)  =  0 
CONTINUE 
CONTINUE 

DO  90  J  =  2,JMAX,2 
DO. 88  I  =  1,IMAX,2 

CRYST3(I,J,K)  =  2 
CONTINUE 
CONTINUE 

DO  100  J  =  1.JMAX.2 
DO  98  I  =  2,IMAX,2 

CRYST3(I,J,K)  =  2 
CONTINUE 
CONTINUE 
CONTINUE 

DO  124  K  =  1,KMAX,2 
DO  114  J  =  1,JMAX,2 
DO  112  I  =  1,IMAX,2 
CRYST3(I,J,K)  =  2 
CONTINUE 
CONTINUE 

DO  123  J  =  2.JMAX.2 
DO  122  I  =  2,IMAX,2 
CRYST3(I,J,K)  =  2 
CONTINUE 
CONTINUE 
CONTINUE 

DO  146  K  =  3,KMAX,2 
DO  134  J  =  2,JMAX,2 
DO  132  I  =  1,IMAX,2 
CRYST3(I,J,K)  =  0 
CONTINUE 
CONTINUE 

DO  144  J  =  1,JMAX,2 
DO  142  I  =  2>IMAX,2 
CRYST3(I,J,K)  =  0 
CONTINUE 
CONTINUE 
CONTINUE 
K  =  1 

DO  158  J  =  2,JMAX,2 
DO  156  I  =  1,IMAX,2 
CRYST3(I,J,K)  =  1 
CONTINUE 
CONTINUE 

DO  172  J  =  1,JMAX,2 
DO  170  I  =  2,IMAX,2 
CRYST3(I,J,K)  =  1 
CONTINUE 
CONTINUE 


END  OF  THREE  DIMENSIONAL  CRYSTAL  SET  UP  AND  THE  BEGINNING  OF  THE  ONE 
DIMENSIONAL  EQUIVELENCE.   THE  ONE  DIMENSIONAL  ARRAY  CRYSTI L  )  VALUES 
ARE  IDENTICAL  TO  THE  THREE  DIMENSIONAL  CRYST3I I , J ,K  ) . 


DO  188  K  =  1,KMAX 
DO  186  J  =  1,JMAX 
DO  184  I  =  1,IMAX 

L  =  I  +( (J-l  1*IMAX)  ♦  ( (K- 
CRYST(L)  =CRYST3(I,J,K) 
CONTINUE 
CONTINUE 
CONTINUE 
KMAX1  =  KMAX  -1 


1)*PLNMAX) 
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c 
c 
c 

c 
c 

00220 
00230 
00240 
00250 
00260 
00270 
00280 
00290 
00300 
00310 
00320 
00330 
00340 
00350 
00360 
00370 
00380 
00390 
00400 
00410 
00420 
00430 
00440 
00450 
00460 
00470 
00480 
00482 
00484 
00490 
00500 
00510 
00520 
00530 
00540 
00550 
00560 
00570 
00580 
00590 
00600 
00610 
00620 
00630 
00640 
00650 
00660 
00670 
00680 
00690 
00700 
00710 
00720 
00730 
00740 
00750 
00760 
00770 
00780 
00790 
00800 
00810 
00820 


THE  FOLLOWING  GROUP  OF  DO  LOOPS 
ARRAY  CRNNBRI L,R).  THIS  ARRAYS 
NEIGHBORS  TO  THE  CRYSTAL  SITE  L 


DEVELOP  THE  TWO  DIMENSIONAL 
VALUES  ARE  THE  TWELVE  NEAREST 


DO  380  K  =  2.KMAX1 
J  =  1 
1  =  1 

L  =  I  ♦ 

CRNNBRI  L,l)  i 
CRNNBRI  L>2)  = 
CRNNBRI L, 3)  ■ 
CRNNBR(L,4)  " 
CRNNBR(L,5)  ■ 
CRNNBRI  L>6)  ■ 
CRNNBRI  L, 7)  i 
CRNNBRI  L>8)  ■■ 
CRNNBRI  L, 9)  ■ 
CRNNBRI  L, 10) 
CRNNBRI  L,. 11) 
CRNNBRI L, 12) 
CONTINUE 

DO  530  K  =  2,KMAX1 
J  =  1 
I  = 


<K-1)*PLNMAX 

IMAX  ♦  K*PLNMAX 

I  ♦  (JMAX-1  l*IMAX  ♦  K*PLNMAX 

1+1  ♦  K*PLNMAX 

I  ♦  IMAX  ♦  K*PLNMAX 

IMAX  ♦  I JMAX-1 )*IMAX  ♦  <K-1)*PLNMAX 

1*1  ♦  I JMAX-1  )*IMAX  ♦  (K-1)*PLNMAX 

1*1  ♦  IMAX  ♦(K-1)*PLNMAX 

IMAX  +  IMAX  ♦  (K-1)*PLNMAX 

IMAX  ♦  <K-2)*PLNMAX 
=  I  ♦  ( JMAX-1  )*IMAX  ♦  (K-2)*PLNMAX 
=  1+1  ♦  (K-2)*PLNMAX 
=  I  ♦  IMAX  ♦  (K-2)*PLNMAX 


IMAX 
L  =  I 


+  <K-1)*PLNMAX 


CONTINUE 
DO  700  K 


CRNNBRI L>1)  = 
CRNNBRI  L,2)  = 
CRNNBRI L, 3)  a 
CRNNBRI  L, 4)  = 
CRNNBRI  L, 5)  = 
CRNNBRI L ,6)  = 
CRNNBRI  L, 7)  = 
CRNNBRI  L, 8)  " 
CRNNBRI  L, 9)  = 
CRNNBRI  L,10) 
CRNNBRI  L, 11) 
CRNNBRI  L, 12) 

=  2.KMAX1 


1-1  +  K*PLNMAX 

I  +  I JMAX-1  )*IMAX  ♦  K*PLNMAX 

1  +  K*PLNMAX 

I  ♦  IMAX  +K*PLNMAX 

1-1  ♦  I JMAX-1  )*IMAX  +  (K-1)*PLNMAX 

1  ♦  I JMAX-1 )*IMAX  ♦  (K-1)*PLNMAX 

1  ♦  IMAX  ♦  <K-1)*PLNMAX 

1-1  ♦  IMAX  ♦  (K-1)*PLNMAX 

1-1  +  <K-2)*PLNMAX 
=  I  ♦  I JMAX-1  )»IMAX  ♦  <K-2)*PLNMAX 
=  1  ♦  <K-2)*PLNMAX 
=  I  ♦  IMAX  ♦  (K-2)*PLNMAX 


J  =  JMAX 


I  = 

L 


IMAX 

=  K*  PLNMAX 
CRNNBRI  L,l)  ■ 
CRNNBRI  L>2)  = 
CRNNBRI L,3)  = 
CRNNBRI  L,4)  = 
CRNNBRI  L, 5)  = 
CRNNBRI L, 6)  = 
CRNNBRI L, 7)  = 
CRNNBRI  L, 8)  = 
CRNNBRI  L, 9)  = 
CRNNBRI  L, 10) 
CRNNBRI  L, 11) 
CRNNBRI  L,  12) 


1-1  ♦  <J-1)*IMAX 
I  ♦  I J-2)*IMAX  ♦ 
1  ♦  I J-l  1*IMAX  ♦ 
I  +  K*PLNMAX 
1-1  +  I J-2  >*IMAX 
1  ♦  I J-2)*IMAX  ♦ 
1  ♦  (K-l  )*PLNMAX 
1-1  ♦  (K-l  )*PLNMAX 
1-1  ♦  I J-l t*IMAX  ♦ 
=  I  +  (J-2)*IMAX  ♦ 
=  1  ♦  (J-l  >*IMAX  ♦ 
=  I  ♦IK-2)*PLNMAX 


♦  K*PLNMAX 
K*PLNMAX 
K*PLNMAX 

♦  (K-l  l*PLNMAX 
(K-l  )*PLNMAX 


(K-2)*PLNMAX 
(K-2  >*PLNMAX 
(K-2)*PLNMAX 


CONTINUE 
DO  850  K 


2,KMAX1 


J  =  JMAX 


I  = 

L 


1 
=  I 


♦  (J-1)*IMAX  +  (K-1)*PLNMAX 


CRNNBRI L,l) 
CRNNBRI L, 2) 
CRNNBRI L, 3) 
CRNNBRI L,4) 
CRNNBRI L, 5) 
CRNNBRI L, 6) 
CRNNBRI L, 7) 
CRNNBRI L,8) 


IMAX  +  I J-l  )*IMAX 
I  ♦  (J-2I*IMAX  ♦ 
1  +  1  ♦  I J-l  )*IMAX 
I  ♦  K*PLNMAX 
IMAX  ♦  I J-2  )*IMAX 
1+1  +  (J-2)*IMAX  ♦ 
1+1  +<K-1)*PLNMAX 
IMAX  +  (K-1)*PLNMAX 


♦  K*PLNMAX 
K#PLNMAX 
♦  K*PLNMAX 


♦  I K-l  )*PLNMAX 
(K-1)*PLNMAX 
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00830 

CRNNBR(L,9)  =  IMAX 

00840 

CRNNBR(L,10)  =  I  ♦ 

00842 

CRNNBR(L,11)  =  1+1 

00844 

CRNNBR(L>12)  =  I  ♦ 

00850 

CONTINUE 

00860 

JMAX1 

=  JMAX  -1 

00870 

IMAX1 

=  IMAX  -1 

00880 

00 

1050  K  =  2>KMAX1 

00890 

1  =  1 

00900 

DO  1040  J  =  2,JMAX1 

00910 

L  =  I  ♦  I J-l  )*IMAX  ♦ 

00920 

CRNNBR(Ul)  =  IMAX 

00930 

CRNNBR(L,2)  =  I  ♦ 

00940 

CRNNBR(L,3)  =  1+1 

00950 

CRNNBR(L,4)  =  I  ♦ 

00960 

CRNNBR(L,5)  =  IMAX 

00970 

CRNNBRIU6)  =  1  +  1 

00980 

CRNNBR(L,7)  =  1*1 

00990 

CRNNBR(L,8)  =  IMAX 

01000 

CRNNBR(L,9)  =  IMAX 

01010 

CRNNBR(l.lO)  =  I  ♦ 

01020 

CRNNBR(L,11)  =  1*1 

01030 

CRNNBR(L,12)  =  I  ♦ 

01040 

CONTINUE 

01050 

CONTINUE 

01060 

DO 

1230  K  =  2,KMAX1 

01070 

J  =  1 

01080 

DO  1220  I  =  2.IMAX1 

01090 

L  =  I  ♦  <K-1)*PLNMAX 

01100 

CRNNBR(L.l)  =  1-1 

OHIO 

CRNNBR(L,2)  =  I  ♦ 

01120 

CRNNBRIU3)  =  1*1 

01130 

CRNNBR(L,4)  =  I  ♦ 

01140 

CRNNBR(L,5)  =  I-.l 

0x150 

CRNNBR(L>6)  =  1*1 

01160 

CRNNBR(L,7)  =  1*1 

01170 

CRNNBR(L,8)  =  1-1 

01180 

CRNNBR(L,9)  =  1-1 

01190 

CRNNBR(L.IO)  ■  I  ♦ 

01200 

CRNNBR(L.ll)  =  1*1 

01210 

CRNNBR(L,12)  =  I  ♦ 

01220 

CONTINUE 

01230 

CONTINUE 

01240 

DO 

1410  K  =  2.KMAX1 

01250 

I  =  IMAX 

01260 

DO  1400  J  =  2,JMAX1 

01270 

L  =  I  ♦  (J-l  )*IMAX  ♦ 

01280 

CRNNBRIL1)  =  1-1 

01290 

CRNNBR(L,2)  =  I  ♦ 

01300 

CRNNBR(L,3)  =  1  ♦ 

01310 

CRNNBR(L,4I  =  I  +< 

01320 

CRNNBRIL.5)  =  1-1 

01330 

CRNNBR(L,6)  =  1  ♦ 

01340 

CRNNBR(L,7)  =  1  ♦  . 

01350 

CRNNBR(L,8)  =  1-1 

01360 

CRNNBR<U9)  =  1-1 

01370 

CRNNBRt  L,10)  =  I  + 

01380 

CRNNBRIL.,11  )  =  1  ♦ 

01390 

CRNNBR(L,12)  =  I  ♦ 

01400 

CONTINUE 

01410 

CONTINUE 

01420 

DO 

1590  K  =  2,  KMAX1 

01430 

J  =  JMAX 

01440 

DO  1580  I  =  2,IMAX1 

01450 

L  =  I  ♦  (J-l  )*IMAX  ♦ 

01460 

CRNNBR(L>1)  =  1-1 

01470 

CRNNBRIL.2)  =  I  + 

01480 

CRNNBR(L,3)  =1+1 

♦  (J-1)*IMAX 
(J-2)*   IMAX 

♦  (J-l >*IMAX 
(K-2)*PLNMAX 


■f  (K-2)*PLNMAX 

♦  (K-2)*PLNMAX 

♦  (K-2)*PLNMAX 


(K- 


1  )*PLNMAX 
(J-l  l*IMAX 


♦  (K*PLNMAX) 


♦ 
(J-2)*IMAX  ♦  K*PLNMAX 
+  (J-1)*IMAX  +  K*PLNMAX 
J*IMAX  +  K*PLNMAX 

♦  (J-2>*IMAX  +  (K-1)*PLNMAX 
♦  (J-2)*IMAX  ♦  (K-1)*PLNMAX 
♦(J*IMAX)  ♦  (K-1)*PLNMAX 

♦  J*IMAX  +  (K-1»*PLNMAX 

♦  (J-1)*IMAX  ♦  <K-2)*PLNMAX 
(J-2)*IMAX  ♦  (K-2)*PLNMAX 
«•  (J-D*IMAX  ♦  (K-2I*PLNMAX 

J*IMAX  ♦  <K-2)*PLNMAX 


♦  K*PLNMAX 
(JMAX-1  )*IMAX  +  K*PLNMAX 

♦  K*PLNMAX 
IMAX  +K*PLNMAX 
+  ( JMAX-1  )*IMAX 

( JMAX-1  )*IMAX 
IMAX  ♦  (K 
IMAX  ♦  (K 


♦  (K-1)*PLNMAX 

♦  <K-1)*PLNMAX 
1  l*PLNMAX 
1  )*PLNMAX 


(K-2)*PLNMAX 
(JMAX-1  )*IMAX  +  1K-2)*PLNMAX 
♦  <K-2)*PLNMAX 
IMAX  +  (K-2  >*PLNMAX 


(K-1)*PLNMAX 
♦  <J-U»IMAX 
(J-2  )*IMAX  ♦ 
(J-l  >*IMAX    + 


+  K»PLNMAX 
(K*PLNMAX) 
K»PLNMAX 


+  (J*IMAX)    ♦    K*PLNMAX 

♦  (J-2)*IMAX  +  <K-1)*PLNMAX 
(J-2)*IMAX  +  (K-l  l*PLNMAX 
J*IMAX  ♦  (K-l  i*PLNHAX 
+(J*IMAXI  +  (K-1I*PLNMAX 

♦  (J-1I*IMAX    ♦    <K-2)*PLNMAX 
(J-2>*IMAX    ♦    (K-2I*PLNMAX 
<J-1)*IMAX    ♦    <K-2)*PLNMAX 
J*IMAX    ♦    (K-2)*PLNMAX 


(K-1I*PLNMAX 

♦  (J-1)*IMAX 
( J-2)*IMAX    + 

♦  (J-1)*IMAX 


♦K»PLNMAX 
K*PLNMAX 
+   K*PLNMAX 
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01490 

CRNNBR(L,4)  =  I  ♦  K*PLNMAX 

01500 

CRNNBR(L,5)  =  1-1  ♦  (J-2 >»IMAX  ♦  <K-1)*PLNMAX 

01510 

CRNNBR(L,6)  =  1*1  ♦  <J-2I*IMAX  ♦  <K-1)*PLNMAX 

01520 

CRNNBR(L,7)  =  1*1  +  (K-1)*PLNMAX 

01520 

CRNNBR(L,8)  =  1-1  ♦  <K-1)*PLNMAX 

01540 

CRNNBR(L,9)  =  1-1  +  (J-l  )*INAX  ♦  (K-2)«PLNMAX 

01550 

CRNNBR(L,10)  =  I  ♦  <J-2)»IMAX  ♦  <K-2)*PLNMAX 

01560 

CRNNBR(L,11)  =  1  +  1  ♦  ( J-l  l*IMAX  +  t  K.-2  )*PLNMAX 

01570 

CRNNBR(L,12)  =  I  ♦  (K-2)*PLNMAX 

01580 

CONTINUE 

01590 

CONTINUE 

01600 

DO 

1780  K 

=  2.KMAX1 

01610 

DO  1770  J  =  2,JMAX1 

01620 

DO  1760  I  =  2.IMAX1 

01630 

L 

=  I  ♦  <J-1)*IMAX  ♦  <K-1>*PLNMAX 

01640 

CRNNBR(L>1)  =  1-1  ♦  <J-1)*IMAX  ♦  K*PLNMAX 

01650 

CRNNBR(L,2)  =  I  ♦  (J-2)*IMAX  ♦  K*PLNMAX 

01660 

CRNNBR(L,3)  =  1*1  ♦  (J-l  )*IMAX  ♦  K*PLNMAX 

01670 

CRNNBR(L,4)  =  I  4  J*IMAX  ♦  K*PLNMAX 

01680 

CRNNBRIL.5)  =  1-1  ♦  <J-2)»IMAX  ♦  (K-1»*PLNMAX 

01690 

CRNNBR(L,6)  =  1*1  +  (J-2)*IMAX  ♦  <K-1>*PLNMAX 

01700 

CRNNBR(L,7)  =  1*1  «•  J*IMAX  ♦  (K-1)*PLNMAX 

01710 

CRNNBR(L,8)  =  1-1  ♦  J*IMAX  ♦  (K-1»*PLNMAX 

01720 

CRNNBR(L,9)  =  1-1  ♦  (J-l  )*IMAX  ♦  (K-2)*PLNMAX 

01730 

CRNNBRIL.10)  =  I  ♦  (J-2I*IMAX  ♦  (K-2)*PLNMAX 

01740 

CRNNBR(L,11)  =  1*1  +  <J-1)*IMAX  ♦  (K-2)»PLNMAX 

01750 

CRNNBR(L,12)  =  I  ♦  J*IMAX  ♦  (K-2)*PLNMAX 

01760 

CONTINUE 

01770 

CONTINUE 

01780 

CONTINUE 

01790 

DO 

1960  J 

=  2.JMAX1 

01800 

K  =  1 

01810 

DO  1950  I  =  2,IMAX1 

01820 

L 

=  I  ♦  <J-1)*IMAX 

01830 

CRNNBR(L»1)  =  1-1  ♦  (J-1)*IMAX  ♦  PLNMAX 

01840 

CRNNBR(L,2)  =  L 

01850 

CRNNBR(L,3>  =  L 

01860 

CRNNBR(L,4)  =  I  ♦  <J-2)*IMAX  ♦  PLNMAX 

01870 

CRNNBR(L,5)  =  L 

01880 

CRNNBR(L,6)  =  L 

01890 

CRNNBR(L,7)  =  I*l  +  ( J-l  )*IMAX  ♦  PLNMAX 

01900 

CRNNBR(L,8)  =  L 

01910 

CRNNBR(L,9)  =  L 

01920 

CRNNBR(L.IO)  =  I  +  J*IMAX  ♦  PLNMAX 

01930 

CRNNBR(L,11)  =  L 

01940 

CRNNBR(L,12)  =  L 

01950 

CONTINUE 

01960 

CONTINUE 

01970 

DO 

2130  I 

=  2,IMAX1 

01980 

K  =  1 

01990 

J  =  I 

02000 

L 

=  I 

02010 

CRNNBR(L.l)  =  1-1  ♦  PLNMAX 

02020 

CRNNBR(L,2)  =  L 

02030 

CRNNBR(L,3)  =  L 

02040 

CRNNBR(L,4)  =  I  ♦  ( JMAX-1  )*IMAX  ♦  PLNMAX 

02050 

CRNNBR(L,5)  =  L 

02060 

CRNNBR(L>6)  =  L 

02070 

CRNNBR(L,7)  =1+1  +PLNMAX 

02080 

CRNNBR(L,8>  =  L 

02090 

CRNNBR(L,9)  =  L 

02100 

CRMNBR(L,10)  =  I  ♦  IMAX  +  PLNMAX 

02110 

CRNNBRI L>11)  =  L 

02120 

CRNNBR(L,12)  =  L 

02130 

CONTINUE 

02140 

DO 

2300  J 

=  2,JMAX1 

02150 

K  =  1 

02160 

I  = 

IMAX 

•  ■" 


02170 

02180 

02190 

02200 

02210 

02220 

02230 

02240 

02250 

02260 

02270 

02280 

02290 

02300 

02310 

02320 

02330 

02340 

02350 

02360 

02370 

02380 

02390 

02400 

02410 

02420 

02430 

02440 

02450 

02460 

02470 

02480 

02490 

02500 

02510 

02520 

02530 

02540 

02550 

02560 

02570 

02580 

02590 

02600 

02610 

02620 

02630 

02640 

02650 

02652 

02654 

02660 

02670 

02680 

02690 

02700 

02710 

02720 

02730  ' 

02740 

02750 

02760 

02770 

02780 

02790 

02800 

02810 

02820 


L  =  I  +  <J-1)*IMAX 
CRNNBRIL,1)  = 
CRNNBR(L>2)  = 
CRNNBRI  L, 3)  = 
CRNNBR(L,4)  = 
CRNNBRI L,5)  = 
CRNNBRI L, 6)  = 
CRNNBRI L,7)  = 
CRNNBRI L, 8)  = 
CRNNBRI L, 9)  = 
CRNNBRI L,10)  = 
CRNNBRI L, 11)  = 
CRNNBRI L, 12)  = 


-1  ♦  <J-1)*IMAX  +PLNMAX 


+  <J-2)*IMAX  +  PLNMAX 


1  +  <J-1)*IMAX  +  PLNMAX 


♦  J*IMAX  +  PLNMAX 


CONTINUE 
DO  2470  I 
K  =  1 
J  = 
L 


=  2,IMAX1 


JMAX 
=  I  ♦ 


(J-1)*IMAX 


CONTINUE 
00  2640  J  •■ 
K    =  1 
1  =  1 

L  i 


CRNNBRI L,l)  = 
CRNNBRI  L,2)  = 
CRNNBRI  L, 3)  = 
CRNNBRI L, 4)  = 
CRNNBRI  L, 5)  = 
CRNNBR(L,6)  = 
CRNNBR(L,7)  = 
CRNNBRI L, 8)  = 
CRNNBRIL.9)  = 
CRNNBRI L, 10)  = 
CRNNBRI  L, 11)  = 
CRNNBRI  L, 12)  = 

=  2,JMAX1 


-1  +  <J-1)*IMAX  ♦  PLNMAX 


♦  <J-2)*IMAX  ♦  PLNMAX 


♦1  ♦  IJ-1)*IMAX  ♦  PLNMAX 


♦  PLNMAX 


:  I  +  IJ-1)*IMAX 
CRNNBRI L,l l=IMAX 
CRNNBRI L, 2)  = 
CRNNBRI L, 3)  = 
CRNNBRI L, 4)  = 
CRNNBRI L, 5)  = 
CRNNBRI L ,6)  = 
CRNNBRI L, 7)  = 
CRNNBRI L, 8)  = 
CRNNBRI L ,9)  = 
CRNNBRI L, 10)  = 
CRNhJBRIL.il)  = 
CRNNBRI  L, 12)  = 


+  (J-1)*IMAX  +  PLNMAX 


+  <J-2)*IMAX  ♦  PLNMAX 


+1  +  (J-1)*IMAX  +  PLNMAX 


♦  J*IMAX  ♦  PLNMAX 


CONTINUE 
1  =  1 
J  =  1 
K  =  1 


CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 
CRNNBRI 


1.1)  = 

1.2)  = 

1.3)  = 

1.4)  = 

1.5)  = 

1.6)  ' 

1.7)  : 

1.8)  = 

1.9)  : 
1,10) 
1,11) 
1,12) 


MAX  ♦  PLNMAX 


♦  <JMAX-1)*IMAX  ♦  PLNMAX 


+1  ♦  PLNMAX 


I  ♦  IMAX  ♦  PLNMAX 

1 

1 


I  =  IMAX 
L  =  IMAX 


CRNNBRI L,l)  = 
CRNNBRI L,2)  = 
CRNNBRI  L, 3)  = 


1-1 

L 

L 


♦  PLNMAX 


02830 

02840 

02850 

02860 

02870 

02880 

02890 

02900 

02910 

02920 

02922 

02924 

02930 

02940 

02950 

02960 

02970 

02980 

02990 

03000 

03010 

03020 

03030 

03040 

03050 

03060 

03070 

03080 

03090 

03100 

03110 

03120 

03130 

03140 

03150 

03160 

03170 

03180 

03190 

C 

C 

C 

C 

C 

C 


I  =  IMAX 
J  =  JMAX 
K  =  1 

L 


1  =  1 

L  > 


1  ♦ 


CRNNBR(L,4) 
CRNNBRI L,5) 
CRNNBRI  L,6  I 
CRNNBRI  L, 7) 
CRNNBRI  L,8) 
CRNNBRI  L, 9) 
CRNNBRI L, 10 
CRNNBRI L, 11 
CRNNBRI L, 12 


PLNMAX 
CRNNBRI L,l > 
CRNNBRI L, 2) 
CRNNBRI L, 3) 
CRNNBRI L, 4) 
CRNNBRI  L, 5) 
CRNNBRI L, 6 ) 
CRNNBRI L, 7) 
CRNNBRI t, 8) 
CRNNBRI L, 9) 
CRNNBRI L, 10 
CRNNBRI L, 11 
CRNNBRI L, 12 

<J-1)*IMAX 
CRNNBRI L,l) 
CRNNBRI L, 2) 
CRNNBRI L, 3) 
CRNNBRI L, 4) 
CRNNBRI L, 5) 
CRNNBRI L, 6) 
CRNNBRI L, 7) 
CRNNBRI  L, 8) 
CRNNBRI L, 9) 
CRNNBRI L, 10 
CRNNBRI L, 11 
CRNNBRI L> 12 


=  I  ♦  IJMAX-1)*IMAX  ♦  PLNMAX 

=  L 

=  L 

=  1  +  PLNMAX 

=  L 

=  L 

)  =  I  ♦  IMAX  ♦  PLNMAX 
)  =  L 
)  =  L 


■1  ♦  (J-1)*IMAX  ♦  PLNMAX 


♦  (J-2)*IMAX  ♦  PLNMAX 


1  -MJ-1)*IMAX  ♦  PLNMAX 


)  =  I  ♦  PLNMAX 
)  =  L 
)  =  L 


=  IMAX  +  <J-1)*IMAX  ♦  PLNMAX 

=  L 

=  L 

=  I  ♦  (J-2)*IMAX  ♦  PLNMAX 

=  L 

=  L 
=1*1  ♦  <J-1)*IMAX  ♦  PLNMAX 

=  L 

=  L 

)  =  I  ♦  PLNMAX 
)  =  L 
)  =  L 


THE 
AN 


03193 


03195 


03194 
03196 
03200 
03210 
03220 
03230 
] 
03240 
03250 
03260 
03270 
03280 


03282 


FOLLOWING  CODE  IS  USED  TO  SEND  THE  DEVELOPED  ARRAYS  TO 
OUTPUT  FILE. 


PRINT  3193,  "THIS  IS  THE  LATICE  SITES  VS  PLANES' 
FORMAT  I '1' ,47X,A34) 

PRINT  3195,  'SITE  *' ,  'PLANE' 
FORMAT  ( '-' ,10X,A8,3X,A5) 
DO  3194  L  =  l.LMAX 
PRINT  3196,  L,  ZPLNI L  ) 
CONTINUE 

FORMAT  I'  '  ,10X,I8,5X,I3) 
PRINT  3210,  '  NEAREST  NEIGHBOR  TABLE  ' 
FORMAT  I  '1'  ,S3X,A24) 

PRINT  3230,^'  ,1,2,3,4,5,6,7,8,9,10,11,12 

FORMAT  I  ,-',<.X,Al,9X,Il,7X,Il,7X,Il,7X,Il,7X,Il,7X,Il,7X,Il, 
7X,Il,7X,Il,faX,I2,6X,I2,6X,I2) 
LMAX1  =  LMAX  -  PLNMAX 
DO  3270  L  =  1,LMAX1 

PRINT  3280,  L,ICRNNBR(L,R),  R=l,12) 
CONTINUE 

FORMAT  I '  ' ,1318) 
DO  3287  K  =  1,KMAX 

PRINT  3284,  'THIS  IS  PLANE  NUMBER',  K 
DO  3282  J  =  1,JMAX 

PRINT  3286,  ( CRYST3I I ,J,K ) ,  I  =  1,IMAX) 
CONTINUE 


7^ 


03284  FORMAT  ( ' 1 ' ,5X,A20 ,3X, 13) 

03286  FORMAT  (  '  '  ,4I<+ ) 

03287  CONTINUE 

03288  WRITE  ( 1,3289  UMAX, IMAX1 ,JMAX,JMAX1 ,KMAX,KMAX1 ,PLNMAX,LMAX,LMAX1, 
1  ISEED1,ISEED2,MTMSTP,T1,AVGTM,PLSTM, OUTPUT 

03289  FORMAT! '  ',918,'  ',2116,'  ',518) 

DO  3290  L  =  1,LMAX 


WRITE  ( 1,3291  )CRYST(L) 

03290 

CONTINUE 

03291 

FORMAT  ( '  ' ,18) 

DO  3300  L  =  1,LMAX1 

DO  3292  R  =  1,12 

WRITE  (1,3301) 

03292 

CONTINUE 

03300 

CONTINUE 

03301 

FORMAT  ( '  ',318) 

DO  3310  L  =  1,LMAX 

NRITE  ( 1,3311  )ZPLN(L) 

03310 

CONTINUE 

03311 

FORMAT  ( '  *>I8) 

03320 

STOP 

END 

CRNN8R(L,R) 
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APPENDIX    C 
Srir7USE'    CODE    LISTING 
The    following    is    the    computer    code    for    the    diffusion 
simulation    named    *  D I F  F  U  S  E ' . 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


'DIFFUSE' 

MARK  R.  POLNASZEK 

JUNE  1986 


THE  CODE  'DIFFUSE'  HAS  BEEN  DEVELOPED  TO  SIMULATE  THE  DIFFUSION  OF 
AN  ATOM  INTERSTITIALLY  THROUGH  A  FACE-CENTERED-CUBIC  CRYSTAL  LATTICE 
THIS  LATTICE  IS  AN  INPUT  ARRAY  THAT  WAS  DEVELOPED  IN  ANOTHER  PROGRAM 
AND  INPUTED  ALONG  WITH  EACH  SITES  NEAREST-NEIGHBOR,  AND  EACH  SITES 
PLANAR  POSITION.    THE  METHOD  THAT  'DIFFUSE'  USES  TO  MOVE  THE 
DIFFUSING  ATOMS  THROUGH  THE  CRYSTAL  ARRAY  IS  BY  A  MONTE  CARLO  METHOD 
THE  SIMULATION  THEN  DEVELOPS  THE  DIFFUSION  COEFFICIENT  BY  FICK'S 
FIRST  AND  SECOND  LAW. 


THE  FOLLOWING  DEFINES  THE  VARIABLES  AS 
DIMENSIONS  THE  VARIABLE  ARRAYS'  MEMORY 


REAL  OR  INTEGER 
ALLOCATION. 


AND 


INTEGER  PLNMAX,  CRYSTI 1300  )  ,CRYST3( 08,08,16  )  ,CRNNBR( 1300,12 ) ,R, 

1  MAVGTM,  ZPLN11300),AVGTM,TMSTP2,   CONS( 1100 ,16  )  ,T1 ,MTMSTP ,TMSTP, 

2  AVGTM1 ,AVGTM2 ,AT1 ,AVGTM3 , ATMSTP .OUTPUT 

REAL**   RDMK 1300  ),RDM2( 1300), ACONSI 210,16),  AFLUXJ 210 ,16 )  , 

1  ADERC1(210,16)  ,ADERC2( 210,16  )  ,ADERT( 210,16  )  ,DIFC01( 210,16  )  , 

2  DIFC02I 210,16),  FLUX! 1100,16  ) ,  PLNMX1 ,CVOL< 1100 ,16  ) , 

3  DERC1(1100,16),DERC2(1100,16),DERT(1100,16),  ACVOH  210,16 ), 
<*   DERDCJ  210,16  ),DERT2(  210,16) 

READ! 5 ,*  )  IMAX ,IMAX1 , JMAX , JMAX1 ,KMAX ,KMAX1 , PLNMAX , LMAX , LMAX1 , 
1  ISEED1 ,ISEED2 ,MTMSTP ,T1 ,AVGTM ,TMSTP2 , OUTPUT 

DO    10   L  =  1,LMAX 


C 
C 
C 
C 
C 
C 


THE 
THE 


FOLLOWING  CODE  READS  IN  THE  REQUIRED  INPUT  AND  ECHOS 
INPUT  INTO  AN  OUTPUT  FILE  FOR  INSPECTON. 


00010 


00015 
00020 


00030 


00035 

G00<+0 


00045 
00050 
00055 
00056 
00060 


00070 


1.LMAX1 

=  1,12 

(5,*)  CRNNBR(L.R) 


READ  (5,*)  CRYST(L) 
CONTINUE 
DO    20   L  = 
DO    15   R 
READ 
CONTINUE 
CONTINUE 
DO    30   L  =  1,LMAX 

READ  (5,*)  ZPLN(L) 
CONTINUE 

IF  (MTMSTP  .EQ.  0)      GOTO  55 
DO     «+0   TMSTP  =   1,  MTMSTP 
DO     35   K  =  1.KMAX1 

READ  (5,*)  FLUX( TMSTP, K) 
CONTINUE 
CONTINUE 

DO     50   TMSTP  =   1, MTMSTP 
DO     M   K=  1.KMAX1 

READ  15,*)  CONSI TMSTP ,K  ) 
CONTINUE 
CONTINUE 
IF  (OUTPUT  .NE.  1)  GOTO  300 

WRITE  (6,60)  'INPUTED  VALUES' 
FORMAT  (  '1'  ,33X,A1<+) 

WRITE  (6,70)  ' IMAX', 'IMAX1',' JMAX' , ' JMAX1 ' , 'KMAX' , 'KMAX1 ' 
' PLNMAX ' , ' LMAX ' , ' LMAX1 ' , ' TIME  AVG ' , ' TMSTP2 ' 
FORMAT  ( '-' ,11A10) 
WRITE  (6,80)  IMAX, IMAX1, JMAX, JMAX1 , KMAX, KMAX1 , PLNMAX, LMAX, 
LMAX1,AVGTM,TMSTP2 
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00080 


00090 


00100 


00110 


00120 


00130 


00140 


00145 


00146 


00150 
00160 

00170 

00175 


00180 
00190 
C 

c 
c 
c 
c 
c 

00300 


FORMAT  (  *  '  ,11110) 

WRITE  16,90)  'RANDOM  GENERATOR  SEEDS' 
FORMAT  (  '-'  ,29X,A22) 
WRITE  (6,0100)  ISEED1,ISEED2 
FORMAT  (  '  '  ,10X,2I16) 

WRITE  (6,0110)  'TIME  CONSTANTS' 
FORMAT  (  '-'  ,33X,A14) 
WRITE  (6,0120)  'START  TIME',  'END  TIME' 
FORMAT  ( '-' ,10X,2A10) 

WRITE  (6,0130)  T1,MTMSTP 
FORMAT  (  '  '  ,10X,2I10) 
WRITE  (6,0140)   'INPUTED  CRYSTAL' 
FORMAT  ( 'l',32X,A15) 
WRITE  (6,145)  'LATICE  ATOM  -  2 ', 'DIFFUSING  ATOM  -  1' 

'VACANCY  -  0' 
FORMAT  ( ,-,,3A25) 

WRITE  (6,146)  'SITE',  'ATOM  TYPE' 
FORMAT  ( '-',10X,2A10> 
DO     150   L  =  1,LMAX 

WRITE  (6,160)  L.,  CRYST(  L) 
CONTINUE 

FORMAT  ( '  '  ,10X,2I10J 
WRITE  (6,170)  'NEAREST  NEIGHBOR  TABLE' 
FORMAT  ( 'l',29X,A22) 
WRITE  (6,0175)  ' L '  ,1 ,2,3 ,4,5,6,7,8,9,10 ,11 ,12 
FORMAT  ( '-' ,A8,12I8) 
DO    180   L  =  1,LMAX1 

WRITE  (6,190)  L,(CRNNBR(L,R),  R=l,12) 
CONTINUE 

FORMAT   ('  ',1318) 


THE  CODE  BELOW  SETS  UP  CONSTANTS  THAT  ARE  USED  REPEATEDLY 
THROUGHOUT  THE  REMAINDER  OF  THE  PROGRAM. 


PLNMX1   =  2/FLOAT(PLNMAX) 

Tl  =  Tl  ♦  MTMSTP 

MTMSTP  =  MTMSTP  ♦  TMSTP2 

AVGTMl-=  AVGTM  -  99 

ATI  =  INT( FLOAT! T1)/AVGTM)+1 

MAVGTM  =  INTl FLOAT! MTMSTP)/ AVGTM) 


C 
C 

C 

c 
c 
c 
c 
c 

c 
c 
c 
c 
c 
c 
c 
c 


THE  MAJOR 
TIMESTEPS 
ARE  COMPLETED 


THE  BELOW  CODE  IS  THE  BEGINNING  OF 
DO  LOOP  COMPLETES  ALL  THE  REQUIRED 
ALL  CALCULATIONS  OF  THE  SIMULATION 
DO  LOOP. 


DO  1050  ATMSTP  =  ATI, MAVGTM 


THE  FOLLOWING  DO  LOOP  COUNTS  ALL  DIFFUSING  ATOM* 


DO  LOOP.   THIS 
FOR  THE  PROGRAM. 
WITHIN  THIS 


PLANE.  THIS  IS  THE 
ALSO  CALCULATES  THE 
CRYSTAL  PLANE. 


CONCENTRATION 
CONCENTRATION 


OF  EACH  PLANE. 
PER  UNIT  VOLUME 


IN  EACH 
THE  DO  LOOP 
OF  EACH 


00301 
C 


DO  1001   TMSTP  =  1, AVGTM 

DO    301   L  =  1,LMAX 
K  =  ZPLN(L) 
IF  (CRYST(L) 
CONSI TMSTP, K ) 
CVOLI TMSTP, K) 
END  IF 
CONTINUE 


EQ.  1)   THEN 

=CONS( TMSTP, K)  ♦  1 

=  C VOL ( TMSTP, K)  ♦  PLNMX1 
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c 

C    THE  FOLLOWING  CONSTANT  IS  DETERMINED  FOR  USE  IN  THE  OUTPUT  FILE. 

C 

C 

AVGTM2  =  <ATMSTP*AVGTM)-<100  -TMSTP ) 
C 
C 

C    THE  FOLLOWING  IF  STATEMENT  SHORTENS  THE  NUMBER  OF  SITES  THAT  THE 
C    SIMULATION  WILL  LOOK  AT  IF  THE  STATEMENT  IS  TRUE,  AND  ASSIGNS  A 
C   CONSTANT  NUMBER  OF  SITES  IF  THE  STATEMENT  IS  FALSE. 
C 
C 

IF  (AVGTM2.LE.  KMAX1 )  THEN 

INMBR  =  AVGTM2*  PLNMAX 

GOTO  310 

END  IF 

INMBR  =  LMAX1 
C 
C 

C    THE  FOLLOWING  SUB-ROUTINE  CALLED  FROM  THE  IMSL  LIBRARY  SETS  UP 
C    TWO  ARRAYS  OF  UNIFORM  RANDOM  NUMBERS. 
C 
C 

00310  CALL  SRND(ISEED1,RDM1, INMBR, 2,0) 

00311  CALL  SRND(ISEED2,RDM2, INMBR, 2,0) 
C 

C 

C   THE  DO  999  DO  LOOP  LOOKS  AT  EACH  SITE  IN  THE  CRYSTAL  ARRAY 

C    TO  DETERMINE  IF  A  DIFFUSING  ATOM  IS  IN  THAT  SITE.   IF  THERE 

C    IS  A  DIFFUSING  ATOM  IN  THAT  SITE  THE  CODE  MOVES  IT  TO  NEW 

C   SITE  IF  IT  CAN. 

C 

C 

00312  DO    999   LRDM  =  1, INMBR 

00313  L  =  INT(INMBR*RDM11LRDM))  ♦  1 

00314  R  =  INT  (12  *  RDM2ILRDM))  +  1" 
LI  =  CRNNBR(L,R) 

C 

C 

C   THE  FOLLOWING  IF  STATEMENT  DETERMINES  IF  THE  SITE  IS  OCCUPIED 

C    BY  A  DIFFUSING  ATOM.   IF  THE  STATEMENT  IS  FALSE  THE  CODE  ALLOWS 

C   THE  DO  LOOP  TO  CONTINUE . 

C 

C 

IF  (CRYST(L)  .NE.  1)  GOTO  999 
C 
C 

C    THE  FOLLOWING  DO  LOOP  IS  USED  TO  SEND  FIRST  PLANE  ATOMS  TO 
C    THEIR  SPECIAL  CODE. 
C 
C 

00316  IF  (ZPLN(L)  .EQ.  1)   GOTO  320 
C 

C 

C   THE  NEXT  DO  LOOP  MOVES  THE  ATOM  TO  ITS  NEW  SITE  . 

C 

C 

00317  IF  (CRYST(Ll)  .EQ.  0)  GOTO  ( 400, 400, 400 ,400 ,500 ,500 ,500, 

1  500,600,600,600,600)  R 

00318  GOTO  999 
C 

C 

C    THE  FOLLOWING  IF  STATEMENT  MOVES  THE  FIRST  PLANE  ATOMS. 

C 

C 

00320     IF  (CRYST(Ll)  .EQ.  0)   THEN 

FLUXt TMSTP, 1)  =FLUX( TMSTP, 1)  ♦  PLNMX1 

CRYST(Ll)  =  CRYST(L) 
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END 
GOTO 


IF 
999 


C 
C 
C 
C 

c 
c 

00400 
00401 
00403 
00404 
00410 

c 
c 
c 
c 
c 
c 

00500 


THE  FOLLOWING  CODE  MOVES  ALL  ATOMS  OTHER  THAN  FIRST 
IF  THE  NEAREST-NEIGHBOR  SELECTED  WAS  IN  THE  FORWARD 


PLANE 
PLANE. 


K  =  ZPLNIL) 

FLUX(TMSTP,K)  =  FLUXt TMSTP ,K  ) 

IF  (ZPLN(Ll)  .EQ.  KMAX)   GOTO 

CRYST(Ll)  =  CRYST(L) 

CRYST(L)  =  0 

GOTO  999 


♦  PLNMX1 
410 


THE  FOLLOWING 
PLANE. 


CODE  MOVES  the  atom  if  it  is  to  stay  in  the  same 


CRYST(Ll)  = 
CRYST(L)  = 
GOTO  999 


CRYST(L) 


C 
C 
C 
C 
C 
C 
00600 


00999 

01001 

C 

C 

C 

C 

C 

C 


THE  FOLLOWING  CODE  IS  USED 
PLANE  BEHIND  THE  DIFFUSING 


IF  THE  ATOM  IS  TO  MOVE  INTO  THE 
ATOMS  INITIAL  POSITION. 


K  =  ZPLN(Ll) 

FLUX! TMSTP ,K  )  =  FLUXt TMSTP, K  )  -  PLNMX1 

CRYST(Ll)  =  CRYST(L) 


CRYST(L) 
CONTINUE 
CONTINUE 


THE  FOLLOWING  DO 
LAST  TIMESTEP. 


LOOP  FINDS  THE  CONCENTRATION  FOR  THE  END  OF  THE 


01002 

C 

C 

C 
C 
C 
C 


DO  1002  L  =  l.LMAX 

K  =  ZPLNI L  ) 

IF  (CRYST(L)  .EQ. 

CONS(AVGTM+l,K)  = 

CVOH  AVGTM+1,K)    = 

END  IF 
CONTINUE 


1  )  THEN 

CONS! AVGTM*1,K)  ♦  1 

CVOL(AVGTM+l,K)  +  PLNMX1 


THE  FOLLOWING 
DETERMINE  THE 


DO 


DO  LOOPS  FIGURE  THE  DERIVATIVES  NEEDED  TO 
DIFFUSION  COEFFICIENTS. 


CVOLI TMSTP, K-l))/2 


01003 
01004 


01005 
01006 


01007 


1004  TMSTP  =  1,AVGTM 
DO  1003  K  =  2.KMAX1 

DERCK  TMSTP, K)  =< CVOL< TMSTP, K*l ) 
CONTINUE 
CONTINUE 

DO  1006  TMSTP  =  1,AVGTM 
DO  1005  K  =  3,KMAX1-1 
DERC2(  TMSTP, K  )  =  (  DERCK  TMSTP, K  +  D-DERCK  TMSTP, K-l )  )/2 
CONTINUE 
CONTINUE 

DO   1008  TMSTP  ■  1,AVGTM 
DO  1007  K  =  1,KMAX1 

DERTt TMSTP   ,K  )  =  CVOLI TMSTP  +  1 ,K )  -  CVOL( TMSTP, K ) 
CONTINUE 


01008 
C 
C 
C 

c 
c 
c 


01010 
01020 
C 

c 
c 

c 
c 
c 


CONTINUE 


THE  FOLLOWING  DO  LOOP  AVERAGES  THE  DETERMINED  VALUES  AS  NEEDED 
TO  DETERMINE  THE  DIFFUSION  COEFFICIENT. 


DO  1020  TMSTP  =  1,AVGTM 

DO  1010  K  =  1,KMAX1 
ACONSt ATMSTP ,K)  =  ACONSI ATMSTP ,K  )  ♦ 
AF LUX ( ATMSTP, K)  =  AFLUXI ATMSTP ,K  )  ♦ 
ADERCK  ATMSTP, K)  =  ADERCK  ATMSTP  ,K  ) 
ADERC2I ATMSTP, K)  =  ADERC2( ATMSTP ,K  ) 
ADERT( ATMSTP ,K)  =  ADERTI ATMSTP ,K  )  ♦ 
ACVOLI ATMSTP, K)  =  ACVOL( ATMSTP, K  )  ♦ 
CONTINUE 

CONTINUE 


FLOAT! CONS! TMSTP ,K  )  J/AVGTM 
FLUX! TMSTP ,K)  /AVGTM 

♦  DERCK  TMSTP, K)  /AVGTM 

♦  DERC2I TMSTP, K)  /AVGTM 
DERTI TMSTP, K)  /AVGTM 

CVOLI TMSTP ,K )/AVGTM 


THE  FOLLOWING  DO  LOOP 
BY  BOTH  THE  FIRST  LAW 


DO 


DETERMINES  THE  DIFFUSION  COEFFICIENT 
IDJCOl)  AND  SECOND  LAW  (DIFC02)  METHOD. 


01030 

C 

C 

C 

C 

C 

c 
c 


1030  K  =  3,KMAX1-1 
DIFCOK  ATMSTP, K)  = 
DIFC021 ATMSTP ,K )  = 
CONTINUE 


-1*  AFLUXI ATMSTP, K (/ADERCK ATMSTP, K) 
ADERTI ATMSTP ,K I/ADERC2I ATMSTP ,K ) 


THE  FOLLOWING  DO  LOOP  DETERMINES  THE  FIRST 
DIFFUSION  COEFFICIENT  AND  BY  FICK'S  SECOND 
THE  TIME  DERIVATIVE  OF  THE  CONCENTRATION. 


DERIVATIVE  OF  THE 
LAW  (NON-LINEAR) 


01031 
C 
C 
C 

C 

c 


DO  1031  K  =  3,KMAXl-2 

DE RDC I ATMSTP, K)  =  DIFCOK ATMSTP ,K*1 )  -  DIFCOK ATMSTP ,K-1 ) 
DERT2I ATMSTP, K)  =  DERDCI  ATMSTP  ,K  )*ADERCK  ATMSTP  ,K  )  + 
DIFCOK  ATMSTP ,K  )*ADERC2( ATMSTP ,K ) 
CONTINUE 


THE  FOLLOWING  DO  LOOP  CLEARS  THE  CONTENTS  OF  THE  NAMED  ARRAYS. 


IF  (ATMSTP  .EQ.  MAVGTM )  GOTO  1050 
DO  1033  TMSTP  =  1,AVGTM*1 
DO  1032  K  =  1,KMAX 

CONS ( TMSTP, K)  =  0 
FLUXI TMSTP, K)  =  0 
DERCK  TMSTP  ,K)  =  0 
DERC2I TMSTP, K)  =  0 


DERTI TMSTP, K) 
CVOLI TMSTP, K) 


=  0 


01032 

01033 

01050 

C 

C 

C 

C 

C 

C 


CONTINUE 
CONTINUE 
CONTINUE 


THE  FOLLOWING  CODE  WRITES  THE  GENERATED  VALUES  OF  THE  SIMULATION 
TO  AN  OUTPUT  FILE  FOR  ANALYSIS. 


01051 
01060 


01070 


GOTO  I  1051, 1191, 1241), OUTPUT 

WRITE  (6,1060)  'OUTPUTED  VALUES' 
FORMAT  I '1' ,33X,A15) 

WRITE  (  6 , 1070 ) ' IMAX ' , ' IMAX1 ' , ' JMAX ' , ' JMAX1 ' 

•PLNMAX' , ' LMAX ' , ' LMAX1 ' 
FORMAT  (  ' - ' , 9A8  )  » 


'KMAX' ,'KMAXl', 


/  / 


01080 

01090 

01100 

OHIO 

01120 

01130 

01140 

] 
01145 

01146 

01150 
01160 

01170 

01175 


01180 
01190 
01191 


01200 
01210 
01220 
01230 
01240 
01241 


1 
2 


1 
2 

3 


01340 
01350 
01360 
01370 
01380 


HRITEI 6,1080 )IMAX ,IMAX1 ,JMAX, JMAX1 ,KMAX,KMAX1 ,PLNMAX,LMAX, 

LMAX1 
FORMAT  ( '  ' ,918) 

HRITEI 6, 1090  )' RANDOM  GENERATOR  SEEDS' 
FORMAT  ( '-' ,29X,A22) 
HRITE  (6,1100)  ISEE01,ISEED2 
FORMAT  (  '  ' ,10X,2I16  ) 

HRITE  (6,1110)  'TIME  CONSTANTS' 
FORMAT  ( '-' ,33X,A14) 
HRITE  (6,1120)  'START  TIME',  'END  TIME' 
FORMAT  ( '-' ,10X,2A10) 

HRITE  (6,1130)  T1,MTMSTP 
FORMAT  ( '  ' ,10X,2I10) 
HRITE  (6,1140)   'INPUTED  CRYSTAL' 
FORMAT  <  'l',32X,A15) 
HRITEI 6,1145)  'LATICE  ATOM  -  2' , 'DIFFUSING  ATOM  -  1', 

'VACANCY  -  0' 
FORMAT  (  '-'  ,3A25) 

HRITE (6,1146)  'SITE',  'ATOM  TYPE' 
FORMAT  ( '-' ,10XX2A10) 
DO    1150   L  =  1,LMAX 

HRITEJ 6,1160)  L,  CRYST(L) 
CONTINUE 

FORMAT  I '  ',10X,2I10) 
HRITEI  6,11,70)  "NEAREST  NEIGHBOR  TABLE' 
FORMAT  I '1' ,29X,A22) 
HRITE  (6,1175)  'L '  ,1,2,3 ,4,5,6,7,8,9,10,11,12 
FORMAT  ( '-',A8,12I8) 
DO   1180   L  =  1.LMAX1 

HRITE(6,1190)  L,(CRNNBR(L,R),  R=l,12) 
CONTINUE 

FORMAT   ('  ',1318) 
DO   1210   TMSTP  =  AVGTM1,AVGTM 

AVGTM3  =  ( MAVGTM*AVGTM  )  -  ( 100-TMSTP ) 

HRITE  (6,1220)   'PARAMETERS  FOR  TIMESTEP  f, 

HRITE   (6,1230)  'PLANE',  'CON   T-','C/VOLT- 

'1ST  DERV','2ND  DERV','TIME  DERV,  "CON 
DO   1200   K  =  1,KMAX1 

HRITE   (6,1240)   K ,CONS( TMSTP ,K  )  ,CVOL( TMSTP ,K ) , 
FLUXI TMSTP ,K ) ,   DERCK  TMSTP ,K  )  ,DERC2( TMSTP ,K  )  ,DERT( TMSTP ,K ) , 

C0NS(TMSTP+1,K),  CVOLl TMSTP  +  1,K ) 
CONTINUE 
CONTINUE 
FORMAT  ( 
FORMAT  ( 
FORMAT  ( 


AVGTM3 
' ,'FLUX' , 
T+','C/VOL  T*' 


,A23,I6) 
,9A10) 

,2I10,5F10.5,I10,F10.5) 
DO  1350  ATMSTP  =  1  ,MAVGTM 

HRITEI 6,1360)  'AVERAGE  PARAMETERS  FOR  TIMESTEPS', 
I IATMSTP-1  )*AVGTM)  +  l,'TO' ,ATMSTP*AVGTM 
HRITE  (6,1370)  ' PLANE ','AVG  CONC'.'AVG  C/VOL','AVG  FLUX', 

'1ST  DER',     '2ND  DER',     'TIME  DER', 'DIF  COEF  1', 
•DIF  COEF  2',  'DIF  CO  DER',  'TIME  DER  ?' 

DO  1340  K  =  1.KMAX1 
HRITE  (6,1380)  K ,ACONS( ATMSTP ,K  )  ,ACVOL I ATMSTP ,K ) , 
AFLUX  I ATMSTP ,K  )  ,ADERC1( ATMSTP ,K  )  ,ADERC2I ATMSTP ,K  )  , 
ADERT  I ATMSTP ,K  )  ,DIFC01 1 ATMSTP ,K  )  ,DIFC02( ATMSTP ,K ) , 
DERDCI ATMSTP ,K  )  ,DERT2( ATMSTP ,K ) 

CONTINUE 
CONTINUE 
FORMAT  I  '1' ,A31,I5,A2,I5) 
FORMAT  I '-' ,11A11) 
FORMAT  ( '  ' ,I11,10F11.5) 
STOP 
END 
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\   PD! 


"1 


C  A  *  CD 


LE  IKPUT  FOR  %I)IFSET' 


The  following  is  the  items  that  need  to  be  inputed  into 
the  code  *  D I F  S  E  T ' .  The  r.  e  t  h  o  d  to  input  this  data  is  to  put 
it  in  the  same  order  as  shown  without  the  titles.  There  is 
no  format  to  the  length  of  each  data,  but  insure  that  there 
are  t  w o  blanks  between  each  input. 
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IMAX 

JMAX 

KMAX 

MTMSTP 

Tl 

ISEED1 

ISEED2 

04 

04 

06 

00 

01 

447586930 

88<*75732 

AVGTM 

OUTPUT 

100 

1 

}fl 


\  -O  7)  T7  ■ 


sai:ple  input  fc::   xrir^s:^' 

The  folio  wing  input  for  the  code  *  D I F  F  U  S  E '  is  created 
from  the  output  file,  FILE  FT01F001,  from  the  code  'DIFSET'. 
This  is  accomplished  by  renaming  the  file  as  I N  P  0  L  2  !  .'■  T .'. . 
The  titles  seen  in  the  sample  are  not  a  part  of  the  actual 
i  n  d  u  t . 


r.  I 


Sample  input  of  the  constants  required  for  'diffuse' . 
The  sample  input  contains  titles  that  would  not  be  in  the 
actual  input  data. 


IMAX     IMAX1    JMAX     JMAX1    KMAX     KMAX1   PLNMAX   LMAX     LMAX1 
4       3       4       3       6       5      16       96      80 


RANDOM  NUMBER  GENERATOR  SEEDS 
44  7586950 
88475732 


MTMSTP    Tl     AVGTM   TMSTP2   OUTPUT 
0       1      100      100       1 
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Sample  input  for  the  CRYST(L)  array. 


2 

1 
2 
1 
1 
2 
1 
2 
2 
1 
2 
0 
2 
0 
0 
2 
0 
2 
2 
0 


2 
0 
2 
0 
2 
2 
0 
2 
0 
0 
2 
0 
2 
2 
0 
2 
0 


Sample  input  for  the  CRNN8R(L,R)  array. 


20 


18 


1; 

2 
2 
30 
Z 
Z 


63 
60 
61 
52 


<>, 


Sanple    input    for   the  ZPLN(L)   array. 


1 
1 


1 
2 

2 


6 

6 


n  ^ 


APPENDIX  F 
SAMPLE  OUTPUT  OF  *DIFSET' 
The  sample  output  shown  in  this  Appendix  is  the  soir.  e  of 
the  output  used  to  check  the  code  for  accuracy.   This  is  not 
all  of  the  output  generated,  but  a  representation  of  what  is 
generated . 
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THIS  IS  THE  LATICE  SITES  VS  PLANES 


SITE  « 
1 
2 
3 
4 
5 
i 
7 
8 
» 
10 
11 
12 
IS 
l< 
15 
If 
17 
18 
1» 
20 
21 
22 
23 
24 
25 
2( 
27 
28 
2» 
JO 
31 
32 
33 
34 
35 
34 
37 
38 
3» 
40 
41 
42 
43 
44 
45 
44 
47 
48 
4» 
50 
51 
52 
53 
54 
55 
54 
57 
58 
5» 
40 
41 
42 
43 
44 


PLANE 


2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 


NEAREST  NEIGHBOR  TABLE 


1 

2 

5 

« 

5 

4 

7 

8 

9 

10 

11 

12 

1 

20 

1 

1 

29 

1 

1 

18 

1 

1 

21 

1 

1 

2 

17 

2 

2 

SO 

2 

2 

19 

2 

2 

22 

2 

2 

J 

It 

3 

3 

51 

3 

3 

20 

3 

3 

25 

3 

5 

4 

19 

« 

4 

52 

4 

4 

17 

4 

4 

24 

4 

4 

5 

24 

5 

s 

17 

5 

5 

22 

5 

S 

25 

S 

S 

4 

21 

i 

4 

18 

4 

4 

23 

4 

4 

26 

4 

6 

; 

22 

7 

7 

19 

7 

7 

24 

7 

7 

27 

7 

7 

8 

23 

8 

8 

20 

8 

8 

21 

8 

8 

28 

8 

8 

♦ 

28 

9 

9 

21 

9 

9 

24 

9 

9 

29 

9 

9 

10 

25 

10 

10 

22 

10 

10 

27 

10 

10 

30 

10 

10 

11 

24 

II 

11 

25 

11 

11 

28 

11 

11 

31 

11 

11 

12 

27 

12 

12 

24 

12 

12 

25 

12 

12 

32 

12 

12 

15 

52 

13 

15 

25 

13 

13 

30 

13 

13 

17 

15 

15 

14 

29 

14 

14  " 

24 

14 

14 

51 

14 

14 

18 

14 

14 

15 

50 

15 

IS 

27 

IS 

15 

52 

IS 

IS 

19 

IS 

15 

It 

51 

14 

14 

28 

14 

14 

29 

14 

14 

20 

14 

14 

17 

54 

45 

54 

57 

32 

50 

22 

24 

4 

13 

2 

5 

18 

55 

44 

35 

58 

29 

51 

25 

21 

1 

14 

3 

4 

1* 

34 

47 

34 

59 

50 

52 

24 

22 

2 

15 

4 

7 

20 

35 

48 

35 

40 

51 

29 

21 

23 

3 

14 

1 

8 

21 

40 

53 

58 

41 

20 

18 

24 

28 

8 

1 

4 

9 

22 

37 

34 

59 

42 

17 

19 

27 

25 

S 

2 

7 

10 

25 

38 

35 

40 

45 

18 

20 

28 

24 

4 

3 

8 

11 

24 

39 

34 

57 

44 

19 

17 

25 

27 

7 

4 

5 

12 

25 

44 

37 

42 

45 

24 

22 

50 

32 

12 

5 

10 

13 

24 

41 

38 

45 

44 

21 

25 

51 

29 

9 

4 

11 

14 

27 

42 

39 

44 

47 

22 

24 

52 

30 

10 

7 

J  > 

IS 

28 

43 

40 

41 

48 

25 

21 

29 

51 

11 

8 

♦ 

14 

29 

48 

41 

44 

55 

28 

24 

18 

20 

14 

9 

14 

1 

50 

45 

42 

47 

54 

25 

27 

19 

17 

15 

10 

15 

2 

31 

44 

43 

48 

35 

24 

28 

20 

18 

14 

11 

14 

3 

32 

47 

44 

45 

34 

27 

25 

17 

19 

15 

12 

13 

4 

55 

52 

41 

50 

55 

48 

44 

58 

40 

20 

29 

18 

21 

34 

49 

42 

51 

54 

45 

47 

39 

57 

17 

30 

19 

22 

35 

50 

43 

52 

55 

44 

48 

40 

58 

18 

31 

20 

23 

34 

51 

44 

49 

54 

47 

45 

37 

59 

19 

32 

17 

24 

37 

54 

49 

54 

57 

34 

54 

42 

44 

24 

17 

22 

25 

38 

55 

50 

55 

58 

35 

55 

43 

41 

21 

18 

25 

24 

39 

54 

51 

54 

59 

54 

34 

44 

42 

^2 

19 

24 

27 

60 

55 

52 

55 

40 

35 

33 

41 

43 

25 

20 

21 

28 

41 

40 

55 

58 

41 

40 

38 

44 

48 

28 

21 

26 

29 

42 

57 

54 

59 

42 

37 

39 

47 

45 

25 

22 

27 

30 

43 

58 

55 

40 

45 

58 

40 

48 

44 

24 

23 

28 

31 

44 

59 

54 

57 

44 

39 

37 

45 

47 

27 

24 

25 

32 

45 

44 

57 

42 

49 

44 

42 

54 

34 

32 

25 

30 

17 

44 

41 

58 

45 

SO 

41 

43 

35 

35 

29 

24 

31 

18 

47 

42 

59 

44 

51 

42 

44 

34 

34 

30 

27 

32 

19 

48 

45 

40 

41 

52 

43 

41 

35 

35 

31 

28 

29 

20 

49 

48 

77 

44 

49 

44 

42 

54 

54 

34 

45 

54 

37 

50 

45 

78 

47 

70 

41 

43 

55 

55 

35 

44 

35 

38 

51 

44 

79 

48 

71 

42 

44 

54 

54 

54 

47 

36 

59 

52 

47 

80 

45 

72 

43 

41 

55 

55 

55 

48 

35 

40 

55 

72 

45 

70 

75 

52 

50 

58 

60 

40 

33 

38 

41 

54 

49 

44 

71 

74 

49 

51 

59 

57 

57 

34 

59 

42 

55 

70 

47 

72 

75 

50 

52 

60 

58 

58 

35 

40 

45 

54 

71 

48 

49 

74 

51 

49 

57 

59 

59 

34 

57 

44 

57 

74 

49 

74 

77 

54 

54 

62 

64 

44 

37 

42 

45 

58 

75 

70 

75 

78 

53 

55 

45 

41 

41 

38 

45 

44 

59 

74 

71 

74 

79 

54 

54 

44 

62 

42 

39 

44 

47 

40 

75 

72 

73 

80 

55 

53 

61 

65 

45 

40 

41 

48 

41 

80 

75 

78 

45 

40 

58 

50 

52 

48 

41 

44 

55 

42 

77 

74 

79 

44 

57 

59 

51 

49 

45 

42 

47 

54 

45 

78 

75 

80 

47 

58 

40 

52 

50 

44 

43 

48 

55 

44 

79 

74 

77 

48 

59 

57 

49 

51 

47 

44 

45 

56 

"1  I 

J 


THIS  IS  PLANE  NUMBER 
2    12    1 
12    12 
2    12    1 
12    12 


s  9 


THIS  IS  PLANE  NUMBER 

2    0  2    0 

0    2  0    2 

2    0  2    0 

0    2  0    2 
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APPENDIX  C 

SAMPLE  OUTPUT  OP  XPIPFUSP' 
The  following  output  is  a  representation  of  the  output 
that  *  D I F  F  U  S  E '  generates  to  use  in  the  analysis  of  the 
diffusion  prcess. 


OUTPUTED  VALUES 


IMAX    1MAX1      JHAX    JMAX1      KMAX    KMAX1   PLNMAX     LHAX    LHAX1 
«        S       «        i  t        5       It       U       80 


RANDOM  GENERATOR  SEEDS 
8S«2i%«2S        18U688SU 


TIME  CONSTANTS 


START  TIME   END  TIME 
1        100 


"2 


IMPUTED  CRYSTAL 


LAT1CE  ATOM  -  2         DIFFUSING  ATOM  -  1  VACANCY 


SITE  ATOM  TYPE 


1 

2 

2 

1 

J 

2 

4 

1 

S 

1 

1 

2 

7 

1 

8 

2 

» 

2 

10 

1 

11 

2 

12 

1 

13 

1 

14 

2 

15 

1 

It 

2 

17 

0 

18 

2 

If 

0 

20 

2 

21 

2 

22 

0 

2J 

2 

24 

0 

25 

0 

2< 

2 

27 

0 

28 

2 

29 

2 

JO 

0 

31 

2 

32 

0 

33 

2 

J4 

0 

35 

2 

J6 

0 

37 

0 

38 

2 

3» 

0 

40 

2 

41 

2 

42 

0 

43 

2 

44 

0 

45 

0 

46 

2 

47 

0 

48 

2 

4» 

0 

50 

2 

51 

0 

52 

2 

53 

2 

54 

0 

55 

2 

56 

0 

57 

0 

58 

2 

S» 

0 

60 

2 

61 

2 

62 

0 

63 

2 

64 

0 

o  " 


NEAREST  NEIGHBOR  TABLE 


L 

1 

2 

S 

4 

5 

t 

7 

8 

♦ 

10 

11 

12 

1 

20 

1 

1 

29 

1 

1 

18 

1 

1 

21 

1 

1 

2 

17 

2 

2 

30 

2 

2 

19 

m 

2 

22 

2 

2 

5 

18 

3 

5 

31 

5 

5 

20 

5 

S 

25 

5 

3 

4 

19 

4 

4 

52 

4 

4 

17 

4 

4 

24 

4 

4 

5 

24 

5 

5 

17 

s 

5 

22 

5 

S 

25 

5 

5 

f 

21 

6 

6 

18 

6 

t 

25 

t 

t 

2t 

t 

6 

7 

22 

7 

7 

19 

7 

7 

24 

7 

7 

27 

7 

7 

8 

23 

8 

8 

20 

8 

8 

21 

8 

8 

28 

8 

8 

* 

28 

9 

♦ 

21 

9 

9 

26 

9 

9 

29 

9 

9 

10 

25 

10 

10 

22 

10 

10 

27 

10 

10 

50 

10 

10 

11 

26 

11 

11 

25 

11 

11 

28 

11 

11 

51 

11 

11 

12 

27 

12 

12 

24 

12 

12 

25 

12 

12 

52 

12 

12 

IS 

32 

13 

IS 

25 

IS 

IS 

50 

15 

IS 

17 

15 

13 

14 

2* 

14 

14 

26 

14 

14 

51 

14 

14 

18 

14 

16 

15 

30 

15 

15 

27 

15 

15 

52 

IS 

15 

19 

15 

15 

If 

31 

14 

14 

28 

It 

It 

29 

It 

It 

20 

It 

14 

17 

36 

65 

54 

57 

52 

50 

22 

24 

4 

IS 

2 

5 

18 

33 

44 

55 

58 

29 

51 

25 

21 

1 

14 

5 

6 

I* 

34 

47 

56 

59 

50 

52 

26 

22 

2 

15 

4 

7 

20 

35 

48 

55 

40 

51 

29 

21 

25 

5 

It 

1 

8 

21 

40 

33 

58 

41 

20 

18 

24 

28 

8 

1 

6 

9 

22 

37 

34 

39 

42 

17 

19 

27 

25 

S 

2 

7 

10 

25 

38 

55 

60 

45 

18 

20 

28 

26 

f 

5 

8 

11 

24 

59 

it 

37 

44 

19 

17 

25 

27 

7 

6 

5 

12 

25 

44 

57 

42 

45 

24 

22 

50 

52 

12 

5 

10 

13 

26 

41 

38 

43 

46 

21 

23 

51 

29 

9 

6 

11 

14 

27 

42 

39 

44 

67 

22 

24 

52 

30 

10 

7 

12 

IS 

28 

43 

40 

41 

68 

25 

21 

29 

31 

11 

8 

9 

It 

29 

48 

61 

46 

55 

28 

24 

18 

20 

It 

9 

14 

1 

30 

45 

62 

47 

54 

25 

27 

19 

17 

15 

10 

15 

2 

31 

46 

43 

48 

35 

24 

28 

20 

18 

14 

11 

16 

3 

32 

47 

44 

45 

36 

27 

25 

17 

19 

15 

12 

15 

4 

33 

52 

41 

50 

53 

68 

46 

38 

40 

20 

29 

18 

21 

54 

49 

62 

51 

54 

45 

67 

59 

37 

17 

50 

19 

22 

35 

50 

63 

52 

55 

46 

48 

40 

38 

18 

51 

20 

23 

It 

51 

66 

49 

56 

47 

65 

57 

34 

19 

52 

17 

24 

37 

56 

49 

54 

57 

36 

34 

42 

66 

24 

17 

22 

25 

38 

53 

50 

55 

58 

33 

35 

45 

6  1 

21 

18 

23 

26 

34 

54 

51 

56 

59 

56 

56 

44 

42 

22 

19 

24 

27 

40 

55 

52 

53 

60 

35 

35 

41 

43 

25 

20 

21 

28 

41 

60 

53 

58 

61 

40 

38 

46 

48 

28 

21 

26 

2' 

42 

57 

54 

59 

62 

37 

59 

67 

45 

25 

22 

27 

30 

43 

58 

55 

60 

65 

38 

40 

48 

46 

26 

25 

28 

31 

44 

5* 

56 

57 

66 

39 

37 

45 

67 

27 

24 

25 

32 

45 

64 

57 

62 

49 

44 

42 

54 

36 

52 

25 

50 

17 

46 

61 

58 

65 

50 

41 

43 

35 

35 

29 

26 

51 

18 

47 

62 

59 

66 

51 

42 

44 

36 

54 

30 

27 

32 

19 

48 

63 

60 

61 

52 

43 

41 

55 

55 

51 

28 

29 

20 

49 

68 

77 

66 

69 

64 

62 

54 

56 

56 

45 

34 

37 

50 

65 

78 

67 

70 

61 

65 

55 

55 

33 

46 

35 

38 

51 

66 

79 

68 

71 

62 

66 

56 

54 

36 

67 

56 

39 

52 

67 

80 

65 

72 

63 

61 

55 

55 

35 

68 

33 

60 

53 

72 

65 

70 

75 

52 

50 

58 

60 

40 

55 

38 

41 

54 

69 

66 

71 

76 

69 

51 

59 

57 

57 

34 

39 

42 

55 

70 

67 

72 

75 

50 

52 

60 

58 

38 

35 

40 

43 

56 

71 

68 

69 

76 

51 

49 

57 

59 

39 

36 

37 

44 

57 

76 

69 

74 

77 

56 

56 

62 

66 

44 

37 

42 

45 

58 

73 

70 

75 

78 

55 

55 

65 

61 

41 

38 

43 

46 

59 

74 

71 

76 

79 

54 

56 

64 

62 

42 

39 

44 

47 

60 

75 

72 

75 

80 

55 

55 

41 

65 

43 

40 

41 

48 

61 

80 

73 

78 

65 

60 

58 

50 

52 

48 

41 

46 

35 

62 

77 

74 

79 

66 

57 

59 

51 

69 

45 

42 

67 

54 

63 

78 

75 

80 

67 

58 

60 

52 

50 

46 

45 

48 

35 

64 

79 

74 

77 

68 

59 

57 

49 

51 

47 

44 

45 

36 
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PARAMETERS  FOR  T1MESTEP    100 


PLANE 

CON 

T- 

C/VOL  T- 

FLUX 

1ST  DERV 

2ND  DERV 

TIME  DERV 

CON 

T» 

C/VOL  T« 

1 

8 

1.00000 

0, 

.00000 

0.00000 

0.00000 

0.00000 

8 

1.00000 

2 

6 

0.75000 

0 

00000 

-0.12500 

0.00000 

0.00000 

6 

0.75000 

5 

6 

0.75000 

0. 

12500 

-0.31250 

-0.12500 

-0.12500 

5 

0.62500 

4 

1 

0.12500 

0. 

00000 

-0.37500 

0.12500 

0.12500 

2 

0.25000 

5 

0 

0.00000 

0. 

00000 

-0.06250 

0.00000 

0.00000 

0 

0.00000 
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AVERAGE  PARAMETERS  FOR  T1MESTEP 


1TO   100 


PLANE 

AVG  CONC 

AVG  C/VOL 

AVG  FLUX 

1ST  DER 

2ND  DER 

1 

7.99999 

1 .00000 

0.08500 

0.00000 

0.00000 

2 

1. 009*7 

0.75125 

0.07750 

-0.2J812 

0.00000 

5 

4. 18997 

0.52375 

0.07125 

-0.21*25 

0.02187 

4 

2.54990 

0.318?J 

0.06875 

-0.17437 

0.02844 

S 

1.39999 

0.17500 

0. 06875 

-0.159J7 

0.00000 

TIME  DER  DIF  COEF  1  DIF  COEF  2  DIF  CO  DER  TIME  DER  7 
0.00000  0.00000  0.00000  0.00000  0.00000 
0.00750  0.00000  0.00000  0.00000  0.00000 
0.00625  0.52948  0.19608  0.39427  -0.0747* 
0.00250  0.39427  0.08791  0.00000  0.00000 
0.00000     0.00000     0.00000     0.00000     0.00000 


96 


LIST    OF 


...  T, . ,  r 


A  s  ic  eland,     D .  R . ,     7h  o    Science    a  a  d    l  n  p  i  n  e  e  r  i  n  p    o  .       !  a  ':  c  r  i  a  1  s  , 
pp.    9  5-103,       Brooks /Cole    Engineering    Division,    1924. 


Bakker, 


II . ,         S  t  o  1  w  i  j  k  ,       i ' .  A  .  ,  van      d  e  r      Hoy,       I , .         a  n  ci 

Zuurendonk,      T . J .  ,      "Computer     Simulation     Tor     Material 
A  oulications",      '.,'uc  1 .      3  e  t  a  1 1 ,    p .  9  G  ,    19  7  6. 


\l  a  c  a  n  c 


Bennett,   C .  K .   and  Alder,   3 . J  .  ,   "Persistence  of 

A  o  t  i  o  n  in  Hard  Sphere  Crystals",   J_j_    P  h  y  s  .  9  h  a  m .  Soli  '  s  , 
Vol  32,  np.  2  111-2122,  1971. 


Jowker,  i-i .  and 
1973. 


ing,  D . A  .  ,  Surface  Science  ,  V  o 1  /  1  ,  p 


Collins,   R.  (Professor,  University  of  Salford)  to  Harrison, 
3.7.  (Professor,  II  aval  Postgraduate  School),  7  u  n  d  ;  n  a  r.  t  a  ] 


Studio s  i  n  P  a y  s 1 c  s ,  Private  Communication,  27  Seat 
1935. 


;3  e  r 


C  r  a  n  k  ,      J  . ,       .;.  h  e       i  a  t  "a  e  n  a  t  i  c  s     of     Diffusion  ,      p 


U  ::  i:  ore. 


Hi 


niversity    Press,     197  5. 


3  a  r  !;  e  n  ,  L  .  S  . ,  "Diffusion,  Mobility  and  Their  Interrelation 
through  free  Energy  in  E  i  n  a  r  y  Metallic  S  y  s  t  e  n  s  !:  , 
Transactions  ,  American  Institute  of  Mining  and  M  e  t  a  1  1.  u  r  p  j  c  a  1 
Engineers ,    V o 1 .     17  5,    p .     134,     1 9  4  S . 

de  7  ruin,  II.  J.  and  Murch,  G.E.,  P  h  i  1  o  s  o  p  a  j.  c  a  1  .a  -  a  z  i  n  c  _, 
Vol     27,     p.     1475,     1973. 


3  i  n  s  t  e  i  n  , 


b  e  r   die   von   der   mole  k ular-kinetisc h e n 


T  heorie  der  V/ a  r  m  e- g  e  f  o  r  d  e  r  t  e  B  e  w  e  g  u  n  g  von  i  n  r  u  h  e  n  c  e  n 
Flussigkeiten  suspendieren  Teilchen",  A  nnalen  der  P  h  v  :■.  i  3  , 
Vol.  17,  p.  549,  19  05. 


7  i  c  k  ,   A  .  ,   "  U  b  or  D  i  f  f  u  s  i  o  n  "  ,   Pop  <a  .   A  n  n  .  ,   V  o  1 . 

1335. 


~>  v  . 


r 1 i  n  n ,  P . A .   and 
p .  5  4,  1 9 6  1  . 


iC.anus, 


P  h ysical  Review ,   V  o 1 . 


_  ~> 


G  h  e  z  ,  3  .  and  Langlois,  '. .' .  3.  ,  "More  on  the  Concentration 
Dependence  of  Pick's  Laws",  A  ncr  Lea n  Journal  of  P 1: y  s  i c  s  ,  To 
b o    pub L i  she d    J u 1 y    1 9  3 G  .  ■ 


07 


G  i  r  i  f  a  1  c  o  ,     L .  A  . ,     .'.Loric      -  j  r  ■  r  a  t  ;  o  n     i  n    C.  rystals  ,     E  1  a  L  s  t!  e  1  1 
Publishing    Com pan y  ,     19 64 . 


Gordon,     A.J.        and  Ford,     R.A., 

??.    1  1 4  -  1  1  5  ,    John    '. .'  i  1  e  y    G 


ihe       C  h  enist's       Co  ■:.  pan  ion  , 
Sons, Inc., 19  72. 


G u  v  ,  A .  G . ,  T  r  a  n  s  a  c  t  i  o n  s ,  Japanese  Institute  of  Metalurgy, 
Vol    19,    p.    483,     19  78. 

Guy,  A.G.,  Cooper,  W.D.  and  Poole  R.L.,  Material  S  c  i  once, 
Vol    3,    p.    103,    1977. 

Harrison,  D .  E . ,  Jr.,  LiHiJi_L£  t.i.£il  of.  Zii-L  ziiSJLl.  .^  I.JL±.£._J1 » 
Unpublished  notes  for  PI!  59  11  at  the  Naval  Postgraduate 
School,    1985. 

Hartley,     G.S.,     "Diffusion    and    Swelling    of    High    Polymers,     I", 

Transactions ,     Faraday    Society,     Vol.     4  2,     p .     6 ,     1946. 

King,    G.U.,     Ind  .    £  n  g .    Ccn.  ,    Vol.    42,     p.     2475,     1951. 


5  irkendall  , 


.  w 


"Diffusion       of       Zinc       in       Alp  i1  a 


r  a  n  s  act  i  o  n  s  ,      American     Institute     of     '.'■.  i  n  i  n  g     and 


v  a  n  c  e  d 


b  r  a  s  s    , 

Metallurgical    Engineers,    Vol.     147,    p .     104,     1942. 

Le      Claire,      A  .  D  .  ,      P.]!_v  s_  _ic_a_  i_     C__h£j2  A  s__t  r__v  j_     An 
Treatise ,    Vol    10,    Academic    Press,    p .    261,     1970. 

Manning,     J . R . ,     Diffusion    Kinetics     for     Atoms     i  n    C  r  v  s t a  1 s , 
pp.     2  -  1 0  ,  D  .    'van    G  o  s  t  r  a  n  d    Company,     Inc.,     1968. 


t  a  n  o  ,       C .  ,      "On      the      Relation      Between      the 


i  i  f  f  i!  s  i  o  n 
Coefficients    and    Concentrations    of     Solid     Metals", Journal     or 

Phvsic     (Japan),    Vol    8,     p.     109,     1955. 


Metropolis,     G  .  ,        Rosenbluth,     A  .  V.  .  ,        Rosen  bluth,     G  .     .  , 

Teller,     A.G.,     and     Teller,     P.,     J^    C  hen.     Phys.,     Vol.     21, 

p  .    1  G  .  J  7  ,     1  V  "j  j  . 

[■lurch,      G.E.,      "Ciiciuical     Diffusion     in     highly     Defective 
Solids",     '-hi  1  o  s  o  p  h  i  c  a  1      Magazine    A_ ,     Vol     41,      pp.      157-163, 

I  ■  30. 

G  u  r  c  !  1  ,  G . G .  ,  "  a  s t  T  onic  Transport  in  Soli  d  s  ,  1 9 G 1  . 

[lurch,     0 .  G .     and     iiowick,     A .  S . ,     "Simulation     in     Crystalline 
Solids'',     Diffusion     in    C  rystalline    Solids,     pp.     57  9-427,     19  Pi  4. 


no  r  n 


I!  n 


■  ■  u  r  c  n  ,      b .  ii .     and 

Correlation 

'-aaazi.no    A.    Vol       2,      >p.    673-67 


Calculation    of    the    Diffusion 

arlc    Methods",    ?;tji  osoaM  ,-  -,  1 
y  /  9  . 


5  1  a  v 


n  ,  A.J 


a  n  a 


Underbill,  ?.;;.,     e  r  i  c  a  n 
-- n y s i c s  Vol  52,  p.  376,  1984.  


Journal,  o 


99 


initial  jis'i  i  i  tio::  LIST 


1  .   Defense  Technical  Information  Center 
Caneron  Station 
Alexandria,  Virginia  2  2  3  0  A - 6 1 4  5 

2.  Library,  Code  0142 
naval  Postgraduate  School 

i  5  o  n  t  e  r  e  y  ,  C  a  1  i  i"  o  r*n  i  a  9394  2-50  0  0 

3.  Department  Chairman,  Code  61 
Department  of  Physics 

llaval  Postgraduate  School 
Monterey,  California  93943-5000 

4.  Dr.  D.E.  ilarrison,  Jr.,  Code  61  "x 
Department  of  Physics 

llaval  Postgraduate  School 
Monterey,  California  93943-5000 

5.  Dr.  Lester  In g her,  Code  55 
Department  of  Operations  Research 
llaval  Postgraduate  School 
Monterey,  California  93943-5000 


'  o  .  Conies 


6 .   L  t .  1 1 . R .  Polnaszek,  U  S  N 
0  49  Condon  Dr. 
Charleston,  South  Carolina  29412 


L 


100 


/ 

180  70 


EUDLEY  KNOX  LIBRAE? 
JTAVAI  ~"  SCHOOL 

RNIA  93943-8002 


219318 


Thesis 

P732 

c.l 


Polnaszek 

Aspects  of  simula- 
ting interstitial  dif- 
fusion in  a  face- 
centered-cubic  lattice. 


Thesis 

P732 

c.l 


Polnaozek 

Aspects  of  simula- 
ting interstitial  dif- 
fusion in  a  face- 
centered-cubic  lattice. 


